*! aoeplacebo: 		Tool to run and visualise a geographic placebo test                       
*! Vers. Beta 1.1:	29/01/2020                                               
*! Author:			Reinhard A. Weisser                                                  
*! 					Nottingham Trent University                               
*! 					reinhard.weisser@ntu.ac.uk
*! Dependencies:	spmap, geoinpoly, geodist
                                             

********************************************************************************

program define aoeplacebo 
version 14													/// developped in 15.2, tested for 16 and 14.2

********************************************************************************
*** SYNTAX																	 ***	
********************************************************************************

syntax anything(everything),	   							/// aoe design
                                                            ///
		AREAdata(string)									/// polygon dataset
		CELLdata(string)									/// incident-grid database
		ESTimation(string)									/// estimation syntax
		METHOD(string)										/// degree or geo
		GRID(string)										/// grid cell ID
		TID(string)											/// time ID
		ETIME(string)										/// active operation indicator
		EVariable(string)									/// active effect variable
		[POSition(string)]									///	radius & step												
		[ROTation(string)]									///	rotation bound & step												
		[SCALE(string)]										/// upper scaling bound & step		
		[REPlications(numlist min = 1 max ==1 integer)]		/// No. of replications for permutation test
		[SEED(string)]										/// seed
		[MULTIarea(string)]									/// variable denoting start and end time of time variant areas
		[COMPLEX(string)]									/// more complex AoE (inclduing a few lags, duration, squared)
		[TEMPlacebo(string)]								/// time-shift placebos		
		[EGRIDmatch(string)]								/// link of grid cells and op area (centroid / vertex)
		[ADDSTAT(string)]									/// further model statistics
		[EVREPort(string)]									/// AoE variable used for report
		[OVREPort(string)]									/// other variables used for report
		[TEMPREPort(string)]								/// time-shift placebos report details
		[PERMREPort(string)]								/// bounds for permutation report
		[CELLPOLY(string)]									/// cell polygon vertices for vertex link
		[MAPDB(string)]										/// DB for backgrond map
		[MAPCO(string)]										/// COORD for backgrond map
		[PROGress(string)]									/// display progress (yes / no / detail)
		[OUTPut(string)]									/// output folder
		[OUTName(string)]									/// name suffix of placebo test
		[*]
		
		/* To do
		GEOSETting()										/// Plot the underlying geographic setting (true AoE, grid cells, maps?)
		*/

********************************************************************************
*** INPUT PREPARATIONS & CHECKS												 ***	
********************************************************************************

local design `=subinstr("`: word 1 of `0''",",","",1)'

* Create output folder structure
if "`output'" != "" {
	local rep_check `=subinstr("`output'"," ","",.)'
	if strpos("`rep_check'",",replace") {
		local outdir `=subinstr("`rep_check'",",replace","",.)'
		if substr("`rep_check'",-7,.)=="replace" {
			!rmdir `outdir' /s /q
			}
		}
	if strpos("`rep_check'",",replace")==0 {
		local outdir `output'
		}
	if strpos("`outdir'",":") {
		mata: st_numscalar("dir_check", direxists("`outdir'"))
		if scalar(dir_check) == 0 mkdir "`outdir'"
		}
	if strpos("`outdir'",":")==0 {
		mata: st_numscalar("dir_check", direxists("`c(pwd)'/`outdir'"))
		if scalar(dir_check) == 0 mkdir "`c(pwd)'/`outdir'"
		local outdir "`c(pwd)'/`outdir'"
		}
	local outvis "`outdir'/aoe_visualisations"		
	mata: st_numscalar("diraoe_check", direxists("`outvis'"))
		if scalar(diraoe_check) == 0 mkdir "`outvis'"
	local outres "`outdir'/aoe_results"	
	mata: st_numscalar("diraoe_check", direxists("`outres'"))
		if scalar(diraoe_check) == 0 mkdir "`outres'"	
	local outpoly "`outdir'/aoe_areas"		
	mata: st_numscalar("diraoe_check", direxists("`outpoly'"))
		if scalar(diraoe_check) == 0 mkdir "`outpoly'"	
	}
if "`output'" == "" {
	local outvis "`c(pwd)'/aoe_visualisations"
	!rmdir "`outvis'" /s /q	
	mkdir "`outvis'"
	local outres "`c(pwd)'/aoe_results"		
	!rmdir "`outres'" /s /q	
	mkdir "`outres'"
	local outpoly "`c(pwd)'/aoe_areas"		
	!rmdir "`outpoly'" /s /q	
	mkdir "`outpoly'"	
	}
	
if "`outname'" != "" local outname _`outname'
	
* Colour adjustments for version 14
if `c(stata_version)' >= 14 & `c(stata_version)' < 15 {
	local c_p1  "`"0 109 44"'"
	local c_p5  "`"65 171 93"'"
	local c_p10 "`"199 233 192"'"
	local c_n10 "`"252 146 114"'"
	local c_n5  "`"239 59 44"'"
	local c_n1  "`"203 24 29"'"
	}
if `c(stata_version)' >= 15 {
	local c_p1  "green%85"
	local c_p5  "green%50"
	local c_p10 "green%20"
	local c_n10 "red%20"
	local c_n5  "red%50"
	local c_n1  "red%85"
	local c_damp80 "%80"
	}	

********************************************************************************
* Error messages

* Programme dependency check
foreach prog in spmap geoinpoly geodist grc1leg kdens savesome {
	cap which `prog'
	if _rc==111 {
		di as error "Programme {opt `prog'} not found. Please install {opt `prog'}. "
		exit 199
		}
	}

* Design specification check	
if !inlist("`design'","diagnostic","permutation") {
	di ""
	di as error "Please specify the AoE placebo design ({opt diagnostic} or {opt permutation})."
	exit 198
	}

* Basic AoE specification check	
if "`position'" == "" & "`rotation'" == "" & "`scale'" == "" {
	di ""
	di as error "It's really hard to run an AoE placebo test without any area placebos..."
	di as error "       ... please specify at least one of the three AoE dimensions:"
	di as error "       {opt position()} {opt rotation()} {opt scale()}"
	di ""
	exit 198	
	}

* Position input checks	
local pos_check `=subinstr("`position'"," ","",.)'
if strpos("`pos_check'",",uniform") | strpos("`pos_check'",",") == 0 {
	local p_distr uniform
	}
if strpos("`pos_check'",",centered") local p_distr centered		
if !inlist("`p_distr'","uniform","centered") {
	di ""
	di as error "Please specify a valid suboption for {opt position()}:"
	di as error "       {opt uniform}  {opt centered}."
	di ""
	exit 198	
	}		
if "`position'" != "" {
	if strpos("`position'",",") local pos_range `=substr("`position'",1,strpos("`position'",",")-1)'
	if strpos("`position'",",") == 0 local pos_range `position'
	local pos_r = `: word 1 of `pos_range''
	cap confirm number `pos_r'
	if _rc == 7 | `pos_r' < 0 local pos_error = 1	
	if "`design'" == "diagnostic" {
		local pos_s = `: word 2 of `pos_range''	
		cap confirm number `pos_s'
		if _rc == 7 | `pos_s' < 0 | `pos_s' > `pos_r' | mod(`pos_r',`pos_s')!=0 local pos_error = 1
		}
	if "`pos_error'" == "1"  | "`pos_range" == ""  {
		di ""
		di as error "Please specify valid position parameters."
		exit 125	
		}	
	}
	
* Rotation input checks	
local rot_check `=subinstr("`rotation'"," ","",.)'
if strpos("`rot_check'",",uniform") | strpos("`rot_check'",",")==0 {
	local r_distr uniform
	}
if strpos("`rot_check'",",centered") local r_distr centered		
if !inlist("`r_distr'","uniform","centered") {
	di ""
	di as error "Please specify a valid suboption for {opt rotation()}:"
	di as error "       {opt uniform}  {opt centered}."
	di ""
	exit 198	
	}		
if "`rotation'" != "" {
	if strpos("`rotation'",",") local rot_range `=substr("`rotation'",1,strpos("`rotation'",",")-1)'
	if strpos("`rotation'",",") == 0 local rot_range `rotation'
	local n_rot = `: word count `rot_range''
	local rot_b = `: word 1 of `rot_range''
	cap confirm number `rot_b'
	if _rc == 7 | !inrange(`rot_b',0,180) local rot_error = 1
	if "`design'" == "diagnostic" & "`n_rot'" == "2" {
		local rot_s = `: word 2 of `rot_range''
		cap confirm number `rot_s'
		if _rc == 7 | `rot_s' < 0 | `rot_s' > `rot_b' local rot_error = 1
		}
	if "`rot_error'" == "1"  | "`rot_range" == "" | ("`n_rot'" == "1" & "`design'" == "diagnostic") {
		di ""
		di as error "Please specify valid rotation parameters."
		exit 125	
		}	
	}	
		
* Scale input checks	
local scale_check `=subinstr("`scale'"," ","",.)'
if strpos("`scale_check'",",uniform") | strpos("`scale_check'",",")==0 {
	local s_distr uniform
	}
if strpos("`scale_check'",",centered") local s_distr centered		
if !inlist("`s_distr'","uniform","centered") {
	di ""
	di as error "Please specify a valid suboption for {opt scale()}:"
	di as error "       {opt uniform}  {opt centered}."
	di ""
	exit 198	
	}		
if "`scale'" != "" {
	if strpos("`scale'",",") local scale_range `=substr("`scale'",1,strpos("`scale'",",")-1)'
	if strpos("`scale'",",") == 0 local scale_range `scale'
	local n_scale = `: word count `scale_range''
	if "`design'" == "diagnostic" & `n_scale' == 2 {
		local scale_ub `: word 1 of `scale_range''
		local scale_s `: word 2 of `scale_range''
		cap confirm number `scale_ub'
		if _rc==7 local scale_error = 1	
		cap confirm number `scale_s'
		if _rc==7 local scale_error = 1	
		if "`scale_error'" != "1" &	(`scale_s' < 0 | `scale_s' > `=`scale_ub'-1') local scale_error = 1
		}
	if "`design'" == "diagnostic" & `n_scale' == 3 {
		local scale_lb `: word 1 of `scale_range''
		local scale_ub `: word 2 of `scale_range''
		local scale_s `: word 3 of `scale_range''
		cap confirm number `scale_lb'
		if _rc==7 local scale_error = 1
		cap confirm number `scale_ub'
		if _rc==7 local scale_error = 1	
		cap confirm number `scale_s'
		if _rc==7 local scale_error = 1	
		if "`scale_error'" != "1" {
			if `scale_s' < 0 | `scale_s' > `=`scale_ub'-1' local scale_error = 1
			if `scale_lb' > 1 | `scale_lb' >= `scale_ub' local scale_error = 1
			}
		}
	if "`design'" == "permutation" & `n_scale' == 2 {
		local scale_lb `: word 1 of `scale_range''
		local scale_ub `: word 2 of `scale_range''
		cap confirm number `scale_lb'
		if _rc==7 local scale_error = 1
		cap confirm number `scale_ub'
		if _rc==7 local scale_error = 1		
		if "`scale_error'" != "1" & (`scale_lb' > 1 | `scale_lb' >= `scale_ub') local scale_error = 1
		}
	if "`scale_error'" != "1" & `scale_ub' < 1 local scale_error = 1		
	if "`scale_error'" == "1"  | "`scale_range" == ""  {
		di ""
		di as error "Please specify valid scale parameters."
		exit 125	
		}	
	}

* Permutation report restriction check	
if "`permreport'" != "" & "`design'" == "permutation" {
	forvalues j = 1/3 {
		local f`j' `=substr("`permreport'",1,strpos("`permreport'",")"))' 
		local permreport `=subinstr("`permreport'","`f`j''","",1)'
		foreach d in p r s {
			if strpos("`f`j''","`d'=(") local `d'_bounds `=subinstr(subinstr("`f`j''","`d'=(","",.),")","",.)'
			}
		}
	if "`p_bounds'" != "" {
		local pbl : word 1 of `p_bounds'
		local pbu : word 2 of `p_bounds'
		cap confirm number `pbl'
		if _rc==7 local pberror = 1
		cap confirm number `pbu'
		if _rc==7 local pberror = 1		
		if !_rc & `pbu' > `pos_r' local pberror = 1
		}
	if "`r_bounds'" != "" {
		local rbl : word 1 of `r_bounds'
		local rbu : word 2 of `r_bounds'
		cap confirm number `rbl'
		if _rc==7 local pberror = 1
		cap confirm number `rbu'
		if _rc==7 local pberror = 1		
		if !_rc & `rbu' > `rot_b' local pberror = 1		
		}	
	if "`s_bounds'" != "" {
		local sbl : word 1 of `s_bounds'
		local sbu : word 2 of `s_bounds'
		cap confirm number `sbl'
		if _rc==7 local pberror = 1
		if !_rc & `sbl' < `scale_lb' local pberror = 1		
		cap confirm number `sbu'
		if _rc==7 local pberror = 1		
		if !_rc & `sbu' > `scale_ub' local pberror = 1		
		}
	if "`pberror'" == "1" {
		di ""
		di as error "Please specify valid bounds for {opt permreport}."
		exit 125	
		}	
	}

* Check for missing polygon coordinate dataset	
if "`egridmatch'" == "vertex" & "`cellpoly'" == "" {
	di ""
	di as error "You selected {opt egridmatch(vertex)}."
	di as error "Please specify a valid grid cell polygon coordinate dataset in {opt cellpoly()}."
	exit 198
	}

* Progress report check	
if !inlist("`progress'","yes","no","detail","") {
	di ""
	di as error "Please specify a valid entry for the {opt progress()} report option:"
	di as error "       {opt yes}  {opt no}  {opt detail}"
	di ""
	exit 198
	}

* Checks for multiple time-variant AoEs
if "`multiarea'" != "" {
local mat_check `=subinstr("`multiarea'"," ","",.)'
if strpos("`mat_check'",",anchored") | strpos("`mat_check'",",") == 0 {
	local fixpoint anchored
	}
if strpos("`mat_check'",",joint") local fixpoint joint		
if strpos("`mat_check'",",split") local fixpoint split	
if !inlist("`fixpoint'","anchored","joint","split") {
	di ""
	di as error "Please specify a valid suboption for {opt multiarea()}:"
	di as error "       {opt anchored}  {opt joint}  {opt split}."
	di ""
	exit 198	
	}	
if strpos("`multiarea'",",") local ma_periods `=substr("`multiarea'",1,strpos("`multiarea'",",")-1)'
if strpos("`multiarea'",",") == 0 local ma_periods `multiarea'	
	qui {
	tempvar aux aux2
	use "`areadata'", clear
	local ea_start: word 1 of `ma_periods'
	local ea_end: word 2 of `ma_periods'	
	local merror = 0
	cap confirm var `ea_start'
	if _rc == 111 local merror = 1
	cap confirm var `ea_end'
	if _rc == 111 local merror = 1	
	if `: word count `=substr("`multiarea'",1,`=strpos("`multiarea'",",")-1')'' != 2  | `merror' == 1 {
		di ""
		di as error "You activated the option for time-variant AoEs." 
		di as error "Please specify the two variables in the {opt areadata} dataset,"
		di as error "which indicate for each subarea the respective start and end data."
		di ""
		exit 198
		}
	gen `aux' = _n
	bysort _X (`ea_start'): gen `aux2' = _n if missing(_X)
	sort `aux'
	replace `aux2' = `aux2'[_n-1] if missing(`aux2')
	drop _ID 
	rename `aux2' _ID
	order _ID, first
	levelsof _ID, local(subareas)
	foreach sa of local subareas {
		sum ea_start if _ID == `sa', meanonly
		local ea_start_`sa' = `r(mean)'
		sum ea_end if _ID == `sa', meanonly
		local ea_end_`sa' = `r(mean)'
		}
	drop `aux'	
	save "`areadata'", replace	
	}
	}
	
if "`egridmatch'" == "" local egridmatch centroid	

********************************************************************************
* Define output dimensions	
qui {	
	use "`celldata'", clear
	
	sum `tid'
	local time_range = r(max) - r(min) +1

	geoinpoly y_cent x_cent using "`areadata'", unique	
		rename _ID in_parea										
		qui replace in_parea=1 if !missing(in_parea)
		qui replace in_parea=0 if missing(in_parea)
		gen in_act_parea=(in_parea==1 & `etime'==1) if !missing(`etime')
		lab var in_act_parea "Centroid of (grid) cell in active operation area"
		
		if ("`complex'" != "") {
			xtset `grid' `tid'
			if strpos("`complex'","TS_`evariable'") {
				bysort `grid' (`tid'): gen TS_in_act_parea = sum(in_act_parea) if in_act_parea==1
					replace TS_in_act_parea = 0 if TS_in_act_parea==. & !missing(in_act_parea)		
					lab var TS_in_act_parea "# of cumulative time periods with active AoE"	
				}
			if strpos("`complex'","TC_`evariable'") {
				bysort `grid' (`tid'): gen TC_in_act_parea = sum(in_act_parea) if in_act_parea==1
				replace TC_in_act_parea = 0 if in_act_parea==0
				tempvar aux
				gen `aux' = TC_in_act_parea
					replace `aux' = L1.`aux' if `aux' == 0 & L1.`aux' != .
					replace `aux' = L1.`aux' if `aux' == TC_in_act_parea & L1.`aux' != .
					replace TC_in_act_parea = TC_in_act_parea - `aux' if TC_in_act_parea != 0
					lab var TC_in_act_parea "# of consecutive time periods with active AoE"	
				}
			local est_p = subinstr(`"`estimation'"',"`evariable'","`complex'",1)
			local est_p = subinstr(`"`est_p'"',"`evariable'","in_act_parea",.)
			}

		if ("`complex'" == "") local est_p = subinstr(`"`estimation'"',"`evariable'","in_act_parea",.)
		`est_p'
		
		local di_addstat_error ""
		foreach j of local addstat {
			cap confirm scalar e(`j')
			if _rc !=0 {
				local di_addstat_error `di_addstat_error' {opt `j'}
				}
			}
		local diag N df_m
		cap confirm scalar e(F)
		if !_rc {
			local diag `diag' df_r F p_F
			local go_t 1
			}	
		cap confirm scalar e(chi2)
		if !_rc {
			local diag `diag' chi2 p_chi2
			local go_z 1
			}
		cap confirm scalar e(r2)
		if !_rc {
			local diag `diag' r2
			cap confirm scalar e(r2_a)
			if !_rc {
				local diag `diag' r2_a	
				}
			}	
		if _rc != 0 {
			cap confirm scalar e(r2_p)
			if !_rc {
				local diag `diag' r2_p
				}
			}
		foreach j of local addstat {
			cap confirm scalar e(`j')
			if !_rc {
				local diag `diag' `j'
				}
			}
		local n_stats: word count `diag'
		tokenize `diag'
		forvalues i = 1 / `n_stats' {
			local diagb `diagb' (``i'')
			}	
		matrix b=e(b)
		matrix V=e(V)
		local indepvars : colnames b
			// Speed
			/*local idv=subinstr(subinstr(subinstr(subinstr("`indepvars'","c.","",.),"i.","",.),"##"," ",.),"#"," ",.) 
			local idv=subinstr(subinstr(subinstr(subinstr("`idv'"," _cons","",1),"TC_`evariable'","",.),"TS_`evariable'","",.),"in_act_parea","",.)
			local n_count: word count `idv'
			forvalues i=1/`n_count' {
				local idv=subinstr(subinstr(regexr(regexr("`idv'","[0-9]+b?n?\.",""),"[0-9]+o\.",""),"TS_","",.),"TC_","",.)
				}
			local idv `e(depvar)' `idv' `grid' `tid' `evariable' `etime'
			local load_vars: list uniq idv 
			// Speed*/
		tokenize `indepvars' //tokenize `e(indepvars)'
		local aoe_bs ""
		local x_bs ""
		*local n_b `e(df_m)' //mod
		local n_b : word count `indepvars' //mod
		if substr("`indepvars'",-6,.)==" _cons" local n_b = `n_b'-1	//mod	
		forvalues i = 1 / `n_b' {
			if strpos("``i''","in_act_parea") local aoe_bs `aoe_bs' ``i''
			if strpos("``i''","in_act_parea")==0 local x_bs `x_bs' ``i''
			}
		foreach p in j r s {
			local aoe_bs_`p'p `=subinstr("`aoe_bs'","in_act_parea","in_act_`p'parea",.)'
			}
		local b_seq ""
		local b_seqb ""
		local se_seq ""
		local se_seqb ""
		local pv_seq ""
		local pv_seqb ""
		local aoe_coef_error ""
		local n_aoe_b: word count `aoe_bs'
		forvalues i = 1 / `n_b' {
			local colbase: word `i' of `indepvars' //`e(indepvars)'	
			tokenize `aoe_bs'
			forvalues j = 1 / `n_aoe_b' {	
				if "``j''" ==  "`colbase'" {
					local aoe_lab_`j' "`colbase'"
					local aoe_lab_`j'_ev `=subinstr("`aoe_lab_`j''","in_act_parea","`evariable'",.)'
					local aoe_coef_error `aoe_coef_error' `aoe_lab_`j'_ev' 
					local di_aoe_coef_error `di_aoe_coef_error' {opt `aoe_lab_`j'_ev'} 
					local b_seq `b_seq' b`j' 
					local b_seqb `b_seqb' (b`j')
					local se_seq `se_seq' se`j' 
					local se_seqb `se_seqb' (se`j')
					local pv_seq `pv_seq' pv`j' 
					local pv_seqb `pv_seqb' (pv`j')	
					}		
				}		
			}
		if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
			local xb_seq ""
			local xb_seqb ""
			local xse_seq ""
			local xse_seqb ""
			local xpv_seq ""
			local xpv_seqb ""
			local x_coef_error ""
			local n_x_b: word count `x_bs'
			local n_ovrep: word count `ovreport'		
			forvalues i = 1 / `n_b' {
				local colbase: word `i' of `indepvars' //`e(indepvars)'	
				tokenize `x_bs'
				forvalues j = 1 / `n_x_b' {	
					if "``j''" ==  "`colbase'" {
						local x_lab_`j'_ev "`colbase'"
						local x_coef_error `x_coef_error' `x_lab_`j'_ev' 
						local di_x_coef_error `di_x_coef_error' {opt `x_lab_`j'_ev'} 
					forvalues k = 1 / `n_ovrep' { 
						if "`: word `j' of `x_bs''" == "`: word `k' of `ovreport''" {
							local xb_seq `xb_seq' xb`j' 
							local xb_seqb `xb_seqb' (xb`j')
							local xse_seq `xse_seq' xse`j' 
							local xse_seqb `xse_seqb' (xse`j')
							local xpv_seq `xpv_seq' xpv`j' 
							local xpv_seqb `xpv_seqb' (xpv`j')	
							}		
						}		
					}
				}	
			}	
		}	
	}	
	
* Checks for temporal placebos	
local R_check `=subinstr("`templacebo'"," ","",.)'
if strpos("`R_check'",",empty") | strpos("`R_check'",",")==0 {
	local seriesmargin empty
	}
if strpos("`R_check'",",zero") local seriesmargin zero		
if strpos("`R_check'",",preserve") local seriesmargin preserve
if !inlist("`seriesmargin'","empty","zero","preserve") {
	di ""
	di as error "Please specify a valid suboption for {opt templacebo()}:"
	di as error "       {opt empty}  {opt zero}  {opt preserve}."
	di ""
	exit 198	
	}	
if "`templacebo'" != "" {
	local ts ts
	local tsb (ts)
	if strpos("`templacebo'",",") local Range `=substr("`templacebo'",1,strpos("`templacebo'",",")-1)'
	if strpos("`templacebo'",",") == 0 local range `templacebo'
	if "`Range" == "" {
		di ""
		di as error "Please specify a valid range for the temporal placebo shift factor."
		exit 125	
		}
	local s_start = `: word 1 of `Range''
	local s_end = `: word 2 of `Range''	
	cap confirm integer number `s_start'
	if _rc == 7 local Rerror = 1
	cap confirm integer number `s_end'
	if _rc == 7 local Rerror = 1
	if ("`Rerror'" != "1" & `s_start' >= `s_end') | "`Rerror'" == "1"  {
		di ""
		di as error "Please specify a valid range for the temporal placebo shift factor."
		exit 125	
		}
	local n_tprange = `s_end' - `s_start' + 1
	if abs(`s_start') >= `time_range' | abs(`s_end') >= `time_range' | `n_tprange' >= `time_range' {
		di ""
		di as error "The maximum time-wise shift factor for the temporal AoE placebo"
		di as error		"is larger than the number of observed time periods (in `tid')."
		exit 125	
		}	
	}
	
if "`tempreport'" != "" {
	local r_check `=subinstr("`tempreport'"," ","",.)'
	if strpos("`r_check'",",line") {
		local tpgtype line
		}
	if strpos("`r_check'",",") & strpos("`r_check'",",line")==0 {
		di ""
		di as error "Invalid option for {opt tempreport}."
		exit 198	
		}		
	if strpos("`tempreport'",",") local range `=substr("`tempreport'",1,strpos("`tempreport'",",")-1)'
	if strpos("`tempreport'",",") == 0 local range `tempreport'
	if "`range" == "" {
		di ""
		di as error "Please specify a valid range for the temporal placebo diagnostic report."
		exit 125	
		}
	local sr_start = `: word 1 of `range''
	local sr_end = `: word 2 of `range''	
	cap confirm integer number `sr_start'
	if _rc == 7 local rerror = 1
	cap confirm integer number `sr_end'
	if _rc == 7 local rerror = 1
	if ("`rerror'" != "1" & `sr_start' >= `sr_end') | "`rerror'" == "1"  {
		di ""
		di as error "Please specify a valid range for the temporal placebo diagnostic report."
		exit 125	
		}		
	}	

* Specification check for complex AoE	
if "`complex'" != "" {
	local check = 0
	forvalues j = 1 / `n_aoe_b' {	
		if "`evreport'" == "`: word `j' of `aoe_coef_error''" local check = 1
		}
	if `check' != 1 {
		di ""
		di as error "You selected COMPLEX AoE."
		di as error "Please specify one of the following AoE estimates to be displayed in the diagnostic report:"
		di as error "       `di_aoe_coef_error'"
		di ""
		exit 198
		}
	}

* Further report checks (additional variables)	
if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
	local check2 = 0
	forvalues j = 1 / `n_x_b' {	
		forvalues k = 1 / `n_ovrep' {
			if "`: word `k' of `ovreport''" == "`: word `j' of `x_coef_error''" local check2 = `check2'+1
			}
		}
	if `check2' != `n_ovrep' {
		di ""
		di as error "You selected further explanatory variables to be included in the output."
		di as error "Please specify one of the following estimates to be included in the diagnostic report:"
		di as error "       `di_x_coef_error'"
		di ""
		exit 198
		}
	}
	
* Seed check
if "`design'" == "permutation" & "`seed'" != "" {
	cap confirm integer number `seed' 
	if _rc == 7 {
		di ""
		di as error "Please specify a valid integer number for the seed."
		exit 125	
		}
	}

* Further report checks (additional statistics)	
if "`di_addstat_error'"	!= "" {
	di ""	
	di as error "The following additional statistics are not availabe for {helpb `e(cmd)'}:"
	di as error "       `di_addstat_error'"	
	di ""
	exit 198
	}

* Basic map data check	
if ("`mapdb'" != "" & "`mapco'" == "") | ("`mapdb'" == "" & "`mapco'" != "") {
	di ""	
	di as error "To obtain placebo maps, please specify valid coordinate and database datasets"
	di ""
	exit 198
	}
	
********************************************************************************
*** START DISPLAY															 ***	
********************************************************************************
	
use "`areadata'", clear
qui count if _X==. & _Y==.
if `r(N)'>1 local pol s

di ""
di ""
di in red "********************************   {opt AOE PLACEBO}   ********************************"
di ""

di as text "{opt Design}:            `design'"
di "{opt Layout}:            `r(N)' polygon`pol' with `=_N-2*`r(N)'' vertices"
di "{opt Generation method}: `method' adjustment"
di "{opt Link procedure}:    `egridmatch' based"
di ""
di "*********************************************************************************"
	
********************************************************************************
*** GEODETIC PREPARATIONS													 ***	
********************************************************************************

qui {	
	*** 1.1 General preparations for the geodetic adjustment ***	
	if ("`method'" == "geodetic") | ("`design'" == "permutation") {	
		
	* Add position adjustment to original area
	use "`areadata'", clear
	
	gen aux=_n
	
	if inlist("`fixpoint'","anchored","split") local centgen "bysort _ID `ea_start' `ea_end' (aux):"
	
	if inlist("`fixpoint'","joint") {
		tempvar dupcheck centcheck
		bysort _ID `ea_start' `ea_end' (aux): gen `dupcheck' = _n
		bysort _ID `dupcheck' (aux): gen `centcheck' = _n
		sort aux
		local split_cond "if `centcheck' == 1"
		} 

	* calculate the centroid
	`centgen' gen a=0.5*(_X*_Y[_n+1]-_X[_n+1]*_Y) `split_cond'
	`centgen' egen A=total(a) `split_cond'
		replace A=abs(A) `split_cond'
	`centgen' gen c_x=(_X+_X[_n+1])*(_X*_Y[_n+1]-_X[_n+1]*_Y) `split_cond'
	`centgen' egen x_cent=total(c_x) `split_cond'
		replace x_cent=abs(1/(6*A)*x_cent) `split_cond'
	`centgen' gen c_y=(_Y+_Y[_n+1])*(_Y*_X[_n+1]-_Y[_n+1]*_X) `split_cond'
	`centgen' egen y_cent=total(c_y) `split_cond'
	replace y_cent=abs(1/(6*A)*y_cent) `split_cond'
	
	if "`fixpoint'" == "anchored" {
		* Set anchor centroid (earliest area)
		sum `ea_start'
		gen anchor = 1 if `ea_start' == `r(min)'
		sum x_cent if anchor == 1, meanonly 
		replace x_cent = `r(mean)'
		sum y_cent if anchor == 1, meanonly 
		replace y_cent = `r(mean)'
		local kp anchor
		}
		
	if inlist("`fixpoint'","joint") {
		sum x_cent, meanonly
		replace x_cent = r(mean) if missing(x_cent)
		sum y_cent, meanonly
		replace y_cent = r(mean) if missing(y_cent)
		}

	keep _ID _X _Y ?_cent aux `kp' `ea_start' `ea_end'		
	
	tempname d2r a b f
	scalar `d2r' = c(pi) / 180
	scalar `a' = 6378137
	scalar `f' = 1/298.257223563 // WGS-84 ellipsiod	
	scalar `b' = `a' - `a' * `f'
	
	local maxiter 1000
	
	* Distance between original centroid and original vertices
	geodist y_cent x_cent _Y _X, gen(disto) 
	
	local lon1 x_cent
	local lon2 _X
	local lat1 y_cent
	local lat2 _Y		
	local generate dist_test
	
	* Get radians between original centroid and original vertices
	tempvar touse
	gen `touse' = 1 `if' `in'
	qui replace `touse' = 0 if mi(`lat1',`lon1',`lat2',`lon2')
	qui replace `touse' = 0 if `touse' == .			

	* Implement Vincenty's (1975) inverse solution (drawing on GEODIST)
	
	* First approximation, equation 13
	tempvar L lambda
	gen double `L' =  `d2r' * (`lon2' - `lon1') if `touse'
	gen double `lambda' = `L' if `touse'			
	
	tempvar U1 U2 sin_U1 sin_U2 cos_U1 cos_U2
	gen double `U1' = atan((1-`f') * tan(`d2r' *`lat1')) if `touse'
	gen double `U2' = atan((1-`f') * tan(`d2r' *`lat2')) if `touse'
	gen double `sin_U1' = sin(`U1') if `touse'
	gen double `sin_U2' = sin(`U2') if `touse'
	gen double `cos_U1' = cos(`U1') if `touse'
	gen double `cos_U2' = cos(`U2') if `touse'
	drop `U1' `U2'

	tempvar cont
	gen `cont' = `touse'
	foreach v in sin_sigma cos_sigma sigma sin_alpha ///
		cos_sq_alpha cos_2sigma_m C lambda_old {
		tempvar `v'
		gen double ``v'' = .
		}
		
	local iter 0
	local more 1
	while `++iter' < `maxiter' & `more' {
		// equation 14
		qui replace `sin_sigma' = sqrt((`cos_U2' * sin(`lambda'))^2 + ///
			(`cos_U1' * `sin_U2' - `sin_U1' * `cos_U2' * cos(`lambda'))^2) ///
			if `cont'
		// mark out co-incident points
		qui replace `cont' = 0 if `sin_sigma' == 0
		// equation 15
		qui replace `cos_sigma' = `sin_U1' * `sin_U2' + `cos_U1' * ///
			`cos_U2' * cos(`lambda') if `cont'
		// equation 16
		qui replace `sigma' = atan2(`sin_sigma',`cos_sigma') if `cont'
		// equation 17
		qui replace `sin_alpha' = `cos_U1' * `cos_U2' * sin(`lambda') / ///
			`sin_sigma' if `cont'
		// use trig identity to obtain cos^2 alpha
		qui replace `cos_sq_alpha' = 1 - `sin_alpha'^2 if `cont'
		// equation 18
		qui replace `cos_2sigma_m' = `cos_sigma' - 2 * `sin_U1' * ///
			`sin_U2' / `cos_sq_alpha' if `cont'
		// adjust if both points are on the equator
		qui replace `cos_2sigma_m' = 0 if `cos_sq_alpha' == 0 & `cont'
		// compute new lambda and compare to previous one
		qui replace `lambda_old' = `lambda' if `cont'
		// equation 10
		qui replace `C' = `f' / 16 * `cos_sq_alpha' * ///
			(4 + `f' * (4 - 3 * `cos_sq_alpha')) if `cont'
		// equation 11
		qui replace `lambda' = `L' + (1 - `C') * `f' * `sin_alpha' * ///
		  (	`sigma' + `C' * `sin_sigma' * (`cos_2sigma_m' + ///
		  `C'*`cos_sigma' * (-1 + 2* `cos_2sigma_m'^2))) if `cont'
		// mark out observations that have converged
		qui replace `cont' = 0 if abs(`lambda'-`lambda_old') <= 1e-12
		// we are done if all observations have converged
		sum `cont', meanonly
		local more = r(max)
		}
	
		tempvar u_sq A B delta_sigma
		gen double `u_sq' = `cos_sq_alpha' * (`a'^2 - `b'^2) / (`b'^2) if `touse'
		// equation 3
		gen double `A' = 1 + `u_sq' / 16384 * (4096 + ///
			`u_sq' * (-768 + `u_sq' * (320 - 175 * `u_sq'))) if `touse'
		// equation 4
		gen double `B' = `u_sq' / 1024 * (256 + ///
			`u_sq' * (-128 + `u_sq' * (74 - 47 * `u_sq'))) if `touse'
		// equation 6
		gen double `delta_sigma' = `B' * `sin_sigma' * (`cos_2sigma_m' + ///
			`B' / 4 * (`cos_sigma' * ///
			(-1 + 2 * `cos_2sigma_m'^2) - ///
			`B' / 6 * `cos_2sigma_m' * (-3 + 4 * `sin_sigma'^2) * ///
			(-3 + 4 * `cos_2sigma_m'^2))) if `touse'
		// equation 19; convert to km and then to miles if requested
		gen double `generate' = `b' * `A' * (`sigma' - `delta_sigma') ///
			 / 1000 `km_to_miles' if `touse'
		// co-incident points were marked out of the iteration loop
		qui replace `generate' = 0 if `sin_sigma' == 0 & `touse'
		
		// use an extended missing value to flag observations that failed to converge
		qui replace `generate' = .a if `cont'
		qui count if `generate' == .a
		if r(N) {
			dis as err "Warning: failed to converge due to near-antipodal points"
			dis as err "Replaced distance(s) with missing value .a"
			dis as err "Number of distance(s) affected = " r(N)
			}
		
		// Radians between original centroid and vertices
		gen double rad_vertex=atan2(`cos_U2'*sin(`lambda'),`cos_U1'*`sin_U2'-`sin_U1'*`cos_U2'*cos(`lambda'))
			replace rad_vertex=mod(rad_vertex+2*c(pi), 2*c(pi))
		
		drop `L' `sin_U1' `sin_U2' `cos_U1' `cos_U2' `A' `B' `u_sq' `delta_sigma' ///
		`lambda' `C' `lambda_old' `sigma' `sin_sigma' `cos_sigma' `sigma' ///
		`sin_alpha' `cos_sq_alpha' `cos_2sigma_m' `C' `lambda_old' `cont' `touse'
		
		drop dist_test
		
		tempfile area_geo
		save `area_geo'
	}
	}

********************************************************************************
*** 1 POSITION PLACEBO TEST													 ***	
********************************************************************************
if "`design'" == "diagnostic" {
qui{ 

	*** 1.1 PLACEBO AREA GENERATION ***

if ("`position'" != "") {
	if ("`method'" == "degree") {
	
		if ("`progress'" != "no") {
			nois di " "	
			*if ("`progress'" != "no") nois di in red "*** POSITION AOE PLACEBOS: {opt Degree} adjustments ***"	
			nois di "{opt Placebo type}:      {opt Position AOE placebos}"
			nois di " "	
			}
		clear
		
		local pos_steps = `pos_r'/`pos_s'

		set obs `=`pos_steps'+1'
		gen xj = _n-1
		expand `=`pos_steps'+1'
		bysort xj: gen yj = _n-1

		gen x_adj = xj * `pos_s'
		gen y_adj = yj * `pos_s'
		expand 4

		bysort x_adj y_adj: gen aux = _n
		replace x_adj = x_adj * (-1) if aux == 2
		replace x_adj = x_adj * (-1) if aux == 3
		replace y_adj = y_adj * (-1) if aux == 3
		replace y_adj = y_adj * (-1) if aux == 4

		gen lvec = sqrt((0+x_adj)^2 + (0+y_adj)^2)

		keep if lvec <= `pos_r'
		
		drop aux lvec
		duplicates drop

		sort xj yj	
		gen p_jitter = _n-1

		sum p_jitter
		local jp_max = r(max)
		
		tempfile pos_placebo
		save `pos_placebo'
		
		* Add position adjustment to original area
		use "`areadata'", clear
		
			* calculate the centroid
			gen aux=_n
			
			if inlist("`fixpoint'","anchored","split") local centgen "bysort _ID  `ea_start' `ea_end' (aux):"

			* calculate the centroid
			`centgen' gen a=0.5*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen A=total(a)
				replace A=abs(A)
			`centgen' gen c_x=(_X+_X[_n+1])*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen x_cent=total(c_x)
				replace x_cent=abs(1/(6*A)*x_cent)
			`centgen' gen c_y=(_Y+_Y[_n+1])*(_Y*_X[_n+1]-_Y[_n+1]*_X)
			`centgen' egen y_cent=total(c_y)
			replace y_cent=abs(1/(6*A)*y_cent)	
			
		if "`fixpoint'" == "anchored" {
			* Set anchor centroid (earliest area)
			sum `ea_start'
			gen anchor = 1 if `ea_start' == `r(min)'
			sum x_cent if anchor == 1, meanonly 
			replace x_cent = `r(mean)'
			sum y_cent if anchor == 1, meanonly 
			replace y_cent = `r(mean)'
			local kp anchor
			}
		
		if ("`progress'" != "no") nois _dots 0, title(Generating position placebo areas (Total: `=`jp_max'+1'))		
		forvalues jp = 0/`jp_max' {
			if ("`progress'" != "no") nois _dots `=`jp'+1' 0	
			preserve
			gen p_jitter = `jp'
			merge m:1 p_jitter using "`pos_placebo'", keepusing(?_adj)
				keep if _merge == 3
				drop _merge
			replace x_cent = x_cent + x_adj	
			replace y_cent = y_cent + y_adj	
			replace _X = _X + x_adj	
			replace _Y = _Y + y_adj	
			
			sort p_jitter _ID aux
	
			tempfile pos_parea_`jp'
			save `pos_parea_`jp''
			*save "C:\Users\tew418\Desktop\aoeplacebo_prog_out\AOE_TEST\AOE_MULTI/pos_parea_`jp'.dta", replace	
			restore
			}
			
		use	`pos_parea_0', clear
		forvalues jp = 1/`jp_max' {
			append using "`pos_parea_`jp''"
			}

		gen orig=1 if x_adj==0 & y_adj==0
			replace orig=0 if orig==.
			
		sort p_jitter _ID aux
		tempfile pos_pareas
		save `pos_pareas'
			save "`outpoly'/pos_pareas`outname'.dta", replace	
		}	
			
********************************************************************************

	if ("`method'" == "geodetic") {
		
		use "`area_geo'", clear
	
		if ("`progress'" != "no") {
			nois di " "	
			*if ("`progress'" != "no") nois di in red "*** POSITION AOE PLACEBOS: {opt Geodetic} adjustments ***"	
			nois di "{opt Placebo type}:      {opt Position AOE placebos}"
			nois di " "	
			}
			
		local pos_r = `pos_r'*100 		// max jitter radius (km)
		local pos_s = `pos_s'*100		// jitter step (km)
		local pos_steps = `pos_r'/`pos_s'
		
		tempname d2r a b f
		scalar `d2r' = c(pi) / 180
		scalar `a' = 6378137
		scalar `f' = 1/298.257223563 // WGS-84 ellipsiod	
		scalar `b' = `a' - `a' * `f'
	
		local maxiter 1000
		* Step 1: Get new coordinate with Delta X
		local lon1 x_cent
		local lat1 y_cent
		
		tempvar touse
		gen `touse' = 1 `if' `in'
		qui replace `touse' = 0 if mi(`lat1',`lon1')
		qui replace `touse' = 0 if `touse' == .		

		if ("`progress'" != "no") nois _dots 0, title(Generating position placebo areas (Total: `=(`pos_steps'*2+1)^2', Relevant: approx. `=round(c(pi)/4*(`pos_steps'*2)^2-1,1)'))	
		local count_pos = 0
			
		* DY=0, DX=/=0
		* Rectangular coordinate boundaries, based on four bearings (0, 90, 180, 270 degrees)
		* bound=2: east, bound=4: west	
		foreach recbound of numlist 2 4 {
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(`d2r'*90*(`recbound'-1))
			gen double `cos_alpha1' = cos(`d2r'*90*(`recbound'-1))
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			* Horizontal steps (by jitter step)
			forvalues hps = 1 / `pos_steps' {
				local count_pos = `count_pos' + 1
				if ("`progress'" != "no") nois _dots `count_pos' 0	

				tempvar cont
				gen `cont' = `touse' in 1	
				foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old {
					tempvar `v'
					gen double ``v'' = .
					}
					
				qui replace `sigma' = (`hps'*`pos_s'*1000)/(`b'*`A')		
				
				local iter 0
				local more 1
				while `++iter' < `maxiter' & `more' {
					qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
					qui replace `sin_sigma' = sin(`sigma') if `cont'					
					qui replace `cos_sigma' = cos(`sigma') if `cont'					
					qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
						`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
					qui replace `sigma_old' = `sigma' if `cont'
					qui replace `sigma' = `hps'*`pos_s'*1000/(`b'*`A') + `delta_sigma' if `cont'	

					qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
					sum `cont', meanonly
					local more = r(max)
					}		
			
				tempvar x lambda C L lambda2
				
				if `recbound' == 2 local hs dxp
				if `recbound' == 4 local hs dxn

				gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
				gen double y_cent_`hs'_`hps'= atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
				gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
				gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
				gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
				gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
				gen double x_cent_`hs'_`hps' = `lambda2'*180/c(pi)
				qui replace y_cent_`hs'_`hps' = y_cent_`hs'_`hps'*180/c(pi)
				
				drop `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
					`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont'	
				}
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' 
			}	
			
		rename y_cent y_cent_dxp_0
		rename x_cent x_cent_dxp_0

		reshape long y_cent x_cent, i(aux) string	
		
		split _j, p("_")
		destring _j3, replace
		rename _j3 xstep
		replace xstep=xstep * (-1) if _j2 == "dxn"
		
		drop _j*
		sort xstep
		
		* Step 2: Get new coordinate with Delta Y
		local lon1 x_cent
		local lat1 y_cent
		
		tempvar touse
		gen `touse' = 1 `if' `in'
		qui replace `touse' = 0 if mi(`lat1',`lon1')
		qui replace `touse' = 0 if `touse' == .		
		
		* DY=0, DX=/=0
		* Rectangular coordinate boundaries, based on four bearings (0, 90, 180, 270 degrees)
		* bound=1: north, bound=3: south	
		foreach recbound of numlist 1 3 {
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(`d2r'*90*(`recbound'-1))
			gen double `cos_alpha1' = cos(`d2r'*90*(`recbound'-1))
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			* Horizontal steps (by jitter step)
			forvalues hps = 1 / `pos_steps' {
				if ("`progress'" != "no") {
					forvalues ct = 1 / `=`pos_steps' * 2 + 1' {
						local count_pos = `count_pos' + 1
						nois _dots `count_pos' 0 
						}
					}
		
				tempvar cont
				gen `cont' = `touse' in 1	
				foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old {
					tempvar `v'
					gen double ``v'' = .
					}
					
				qui replace `sigma' = (`hps'*`pos_s'*1000)/(`b'*`A')		
				
				local iter 0
				local more 1
				while `++iter' < `maxiter' & `more' {
					qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
					qui replace `sin_sigma' = sin(`sigma') if `cont'					
					qui replace `cos_sigma' = cos(`sigma') if `cont'					
					qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
						`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
					qui replace `sigma_old' = `sigma' if `cont'
					qui replace `sigma' = `hps'*`pos_s'*1000/(`b'*`A') + `delta_sigma' if `cont'	

					qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
					sum `cont', meanonly
					local more = r(max)
					}		
			
				tempvar x lambda C L lambda2
				
				if `recbound' == 1 local hs dyp
				if `recbound' == 3 local hs dyn

				gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
				gen double y_cent_`hs'_`hps'= atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
				gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
				gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
				gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
				gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
				gen double x_cent_`hs'_`hps' = `lambda2'*180/c(pi)
				qui replace y_cent_`hs'_`hps' = y_cent_`hs'_`hps'*180/c(pi)
				
				drop `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
					`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont'	
				}
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' 
			}			
			
		gen aux2=_n, after(aux)
		
		rename y_cent y_cent_dyp_0
		rename x_cent x_cent_dyp_0

		reshape long y_cent x_cent, i(aux2) string	
		
		split _j, p("_")
		destring _j3, replace
		rename _j3 ystep
		replace ystep=ystep * (-1) if _j2 == "dyn"
		
		drop _j*
		sort xstep ystep
		
		gen x_adj = xstep * `pos_s' 
			lab var x_adj "longitude adjustment (in km)"
		gen y_adj = ystep * `pos_s'
			lab var y_adj "latitude adjustment (in km)"
	
		* Keep centroids in pre-defined radius
		
		sum y_cent if xstep == 0 &  ystep == 0 
		gen y_cent_orig = r(mean)
		sum x_cent if xstep == 0 &  ystep == 0 
		gen x_cent_orig = r(mean)	
		
		geodist y_cent x_cent y_cent_orig x_cent_orig, gen(dist)
		if "`fixpoint'" == "split" {
			bysort xstep ystep: egen mdist = min(dist)
			drop dist
			rename mdist dist
			}	
		keep if dist <= `pos_r'
		
		sort xstep ystep aux
		
		keep _ID ?_cent ?step ?_adj aux disto rad_vertex `kp'		
		
		* Shift the vertices
		rename y_cent y_cent_jgeo
		rename x_cent x_cent_jgeo
		
			local lon1 x_cent_jgeo
			local lat1 y_cent_jgeo

			// at least one lat/lon is a variable; others can be scalar or number
			tempvar touse
			gen `touse' = 1 `if' `in'
			qui replace `touse' = 0 if mi(`lat1',`lon1')
			qui replace `touse' = 0 if `touse' == .					
			
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(rad_vertex)
			gen double `cos_alpha1' = cos(rad_vertex)
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			// iterations
			tempvar cont
			gen `cont' = `touse'	
			foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
				tempvar `v'
				gen double ``v'' = .
				}
			qui replace `sigma' = disto*1000/(`b'*`A')									
			
			local iter 0
			local more 1
			while `++iter' < `maxiter' & `more' {
				qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
				qui replace `sin_sigma' = sin(`sigma') if `cont'					
				qui replace `cos_sigma' = cos(`sigma') if `cont'					
				qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
					`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
				qui replace `sigma_old' = `sigma' if `cont'
				qui replace `sigma' = disto*1000/(`b'*`A') + `delta_sigma' if `cont'	
	
				// mark out observations that have converged
				qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
				// we are done if all observations have converged
				sum `cont', meanonly
				local more = r(max)
				}		
			
			tempvar x lambda C L lambda2

			gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
			gen double _Y = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
			gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
			gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
			gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
			gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
			gen double _X = `lambda2'*180/c(pi)
			qui replace _Y = _Y*180/c(pi)
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
				`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont' `touse'
			cap drop __0000*

		gen a_sort = _n //tmp
		
		gen p_jitter = 0 if xstep == 0 & ystep ==0
		bysort aux p_jitter (a_sort): gen auxj = _n if aux ==1 & p_jitter != 0
		bysort xstep ystep (aux): egen pj = max(auxj)
		replace p_jitter = pj if missing(p_jitter)
	
		drop auxj pj disto
		
		sort p_jitter aux
		
		gen orig=1 if xstep == 0 & ystep == 0
			replace orig=0 if orig==.

		sum p_jitter
		local jp_min = r(min)
		local jp_max = r(max)
		
		tempfile pos_pareas
		save `pos_pareas'
			save "`outpoly'/pos_pareas`outname'.dta", replace	
			
		forvalues jp = `jp_min'/`jp_max' {
			tempfile pos_parea_`jp'
			savesome _all using	`pos_parea_`jp'' if p_jitter == `jp'
			}
		}	// END geo

********************************************************************************
	
	* Identify incidents in various areas
	if ("`progress'" != "no") {
		nois di " "					
		nois di " "	
		nois _dots 0, title(Identifying cells in position placebo areas via `egridmatch' link  (Total: `=`jp_max'+1'))	
		}
	
	if "`egridmatch'" == "centroid" {	
		use "`celldata'", clear
		forvalues jp = 0/`jp_max' {
		if ("`progress'" != "no") nois _dots `=`jp'+1' 0
				if "`multiarea'" == "" {
					geoinpoly y_cent x_cent using "`pos_parea_`jp''", unique 	
					rename _ID in_jparea_`jp'										
					qui replace in_jparea_`jp'=1 if !missing(in_jparea_`jp')
					qui replace in_jparea_`jp'=0 if missing(in_jparea_`jp')
					gen in_act_jparea_`jp'=(in_jparea_`jp'==1 & `etime'==1) if !missing(`etime')
					}
				if "`multiarea'" != "" {
					geoinpoly y_cent x_cent using "`pos_parea_`jp''"	
					gen ea_start = .
					gen ea_end = .
					foreach sa of local subareas {
						replace ea_start = `ea_start_`sa'' if _ID == `sa'
						replace ea_end = `ea_end_`sa'' if _ID == `sa'
						}
					replace _ID = 1 if inrange(`tid',ea_start,ea_end) & !missing(_ID)
					replace _ID = 0 if !inrange(`tid',ea_start,ea_end) | missing(_ID)
					bysort `grid' `tid': egen in_act_jparea_`jp' = max(_ID)
					drop ea_* _ID
					duplicates drop `grid' `tid' in_act_jparea_`jp', force
					}
				lab var in_act_jparea_`jp' "Centroid of (grid) cell in active operation area"
			}
		}
			
	if "`egridmatch'" == "vertex" {		
		if "`multiarea'" == "" {
			use "`cellpoly'", clear
			forvalues jp = 0/`jp_max' {
				if ("`progress'" != "no") nois _dots `=`jp'+1' 0
						
					geoinpoly _Y _X using "`pos_parea_`jp''", unique
					rename _ID aux	
					qui replace aux=1 if !missing(aux)
					qui replace aux=0 if missing(aux)
					bysort `grid': egen in_jparea_`jp'=max(aux)
					drop aux
					}
				
			drop _X _Y
			duplicates drop

			merge 1:m `grid' using "`celldata'"
				keep if _merge==3
				drop _merge
			
			forvalues jp = 0/`jp_max' {
				gen in_act_jparea_`jp'=(in_jparea_`jp'==1 & `etime'==1) if !missing(`etime')
				lab var in_act_jparea_`jp' "Vertex of (grid) cell in active operation area"
				}
			}	
		
		if "`multiarea'" != "" {
			forvalues jp = 0/`jp_max' {
				use "`cellpoly'", clear
				if ("`progress'" != "no") nois _dots `=`jp'+1' 0
				geoinpoly _Y _X using "`pos_parea_`jp''"
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
	
				rename _ID aux	
				qui replace aux=0 if missing(aux)
		
				drop if missing(_X)
				drop _X _Y
				duplicates drop
		
				reshape wide ea_start ea_end, i(`grid') j(aux)
				
				merge 1:m `grid' using "`celldata'", keepusing(`tid')
	
				gen in_act_jparea_`jp' = 0
					lab var in_act_jparea_`jp' "Vertex of (grid) cell in active operation area"

				foreach sa of local subareas {
					cap confirm var ea_start`sa' ea_end`sa'
					if !_rc {
						replace	in_act_jparea_`jp' = 1 if inrange(`tid',ea_start`sa',ea_end`sa') ///
						& !missing(ea_start`sa') & !missing(ea_end`sa')	
						}
					}	
				drop ea_start* ea_end* _merge
				tempfile mat_`jp'
				save `mat_`jp''
				}	
			use "`celldata'", clear
			forvalues jp = 0/`jp_max' {
				merge 1:1 `grid' `tid' using `mat_`jp''
					drop _merge
				}
			}
		}	
		
	tempfile pos_pest	
	save `pos_pest'	
		save "`outres'/pos_pest`outname'.dta", replace
		
	* Optional / Test mapping
	* NOTE: polygon & data can't handle temp files (messed modification)
	if ("`mapdb'" != "") & ("`mapco'" != "") {
		use "`mapdb'", clear
		spmap using "`mapco'", ///
			id(_ID) fcolor(gs12) ocolor(gs12) ///
			polygon(data("`outpoly'/pos_pareas`outname'.dta") by(orig) osize(vthin) ocolor(black cranberry)) ///
			point(data("`outpoly'/pos_pareas`outname'.dta") by(orig) x(x_cent) y(y_cent) ///
			size(vsmall) fcolor(black cranberry))
			graph save "`outvis'/pos_pareas`outname'.gph", replace	
			graph export "`outvis'/pos_pareas`outname'.png", replace	
		}	
	
********************************************************************************

	*** 1.2 PLACEBO AOE ESTIMATIONS ***

	use "`pos_pareas'", clear

		sum p_jitter
		local jp_min = r(min)
		local jp_max = r(max)
		
		tempname memhold
		tempfile pos_pout
		cap postclose `memhold'
		postfile `memhold' p_jitter `ts' `b_seq' `se_seq' `pv_seq' `xb_seq' `xse_seq' `xpv_seq' ///
			`diag' str30(cmd method egridmatch) using "`pos_pout'"
		
	use "`pos_pest'", clear
	
	if ("`progress'" != "no") nois di " "					
	if ("`progress'" != "no") nois di " "					
	if "`templacebo'" != "" {
		if ("`progress'" != "no") nois _dots 0, title(Estimation for position placebo areas ///
		(Total: `=`jp_max'+1' areas x `n_tprange' temp. shifts = `=(`jp_max'+1)*`n_tprange''))
		}
	if "`templacebo'" == ""	{
		if ("`progress'" != "no") nois _dots 0, title(Estimation for position placebo areas (Total: `=`jp_max'+1'))
		}
		
	local s_dot = 0
	forvalues jp = `jp_min'/`jp_max' {
		preserve

		if "`templacebo'" != "" {
			xtset `grid' `tid'

			tempvar aux_treated aux_t1 aux_tT t1 tT

			bysort `grid': egen `aux_treated' = max(in_act_jparea_`jp')

			sum `tid'
			local tmin = r(min)
			local tmax = r(max)
			gen `aux_t1' = in_act_jparea_`jp' if `tid' == `tmin'
			gen `aux_tT' = in_act_jparea_`jp' if `tid' == `tmax'
			bysort `grid': egen `t1' = max(`aux_t1')
			bysort `grid': egen `tT' = max(`aux_tT')
			
			local shift_suffix ""
			forvalues s = `s_start' / `s_end' {
				if `s' < 0 {
					gen in_act_jparea_`jp'_sn_`=`s' * (-1)' = F`=`s' * (-1)'.in_act_jparea_`jp'
					local shift_suffix `shift_suffix' _sn_`=`s' * (-1)'
					if "`seriesmargin'" == "zero" & !missing(in_act_jparea_`jp') ///
						replace in_act_jparea_`jp'_sn_`=`s' * (-1)' = 0 if missing(in_act_jparea_`jp'_sn_`=`s' * (-1)')
					if "`seriesmargin'" == "preserve" & !missing(in_act_jparea_`jp') {
						replace in_act_jparea_`jp'_sn_`=`s' * (-1)' = 1 if missing(in_act_jparea_`jp'_sn_`=`s' * (-1)') & `aux_treated' == 1 & `tT' == 1
						replace in_act_jparea_`jp'_sn_`=`s' * (-1)' = 0 if missing(in_act_jparea_`jp'_sn_`=`s' * (-1)') & (`aux_treated' == 0 | `tT' == 0)
						}
					}
				if `s' > 0 {
					gen in_act_jparea_`jp'_sp_`=`s'' = L`s'.in_act_jparea_`jp'
					local shift_suffix `shift_suffix' _sp_`=`s''
					if "`seriesmargin'" == "zero" & !missing(in_act_jparea_`jp') replace in_act_jparea_`jp'_sp_`=`s'' = 0 if missing(in_act_jparea_`jp'_sp_`=`s'')
					if "`seriesmargin'" == "preserve" & !missing(in_act_jparea_`jp') {
						replace in_act_jparea_`jp'_sp_`=`s'' = 1 if missing(in_act_jparea_`jp'_sp_`=`s'') & `aux_treated' == 1 & `t1' == 1
						replace in_act_jparea_`jp'_sp_`=`s'' = 0 if missing(in_act_jparea_`jp'_sp_`=`s'') & (`aux_treated' == 0 | `t1' == 0)
						}
					}
				}
			}	
			
		foreach shift in "" `shift_suffix' {
			local s_dot = `s_dot' + 1
			if "`shift'" == "" scalar ts = 0
			if strpos("`shift'","_sn_") scalar ts = `=subinstr("`shift'","_sn_","-",1)'
			if strpos("`shift'","_sp_") scalar ts = `=subinstr("`shift'","_sp_","",1)'
			if ("`complex'" != "") {
				xtset `grid' `tid'
				
				if strpos("`complex'","TS_`evariable'") {
					bysort `grid' (`tid'): gen TS_in_act_jparea_`jp'`shift' = sum(in_act_jparea_`jp'`shift') if in_act_jparea_`jp'`shift'==1
						replace TS_in_act_jparea_`jp'`shift' = 0 if TS_in_act_jparea_`jp'`shift'==. & !missing(in_act_jparea_`jp'`shift')		
						lab var TS_in_act_jparea_`jp'`shift' "# of cumulative time periods with active AoE"	
					}
				if strpos("`complex'","TC_`evariable'") {
					bysort `grid' (`tid'): gen TC_in_act_jparea_`jp'`shift' = sum(in_act_jparea_`jp'`shift') if in_act_jparea_`jp'`shift'==1
					replace TC_in_act_jparea_`jp'`shift' = 0 if in_act_jparea_`jp'`shift' == 0
					tempvar aux
					gen `aux' = TC_in_act_jparea_`jp'`shift'
						replace `aux' = L1.`aux' if `aux' == 0 & L1.`aux' != .
						replace `aux' = L1.`aux' if `aux' == TC_in_act_jparea_`jp'`shift' & L1.`aux' != .
						replace TC_in_act_jparea_`jp'`shift' = TC_in_act_jparea_`jp'`shift' - `aux' if TC_in_act_jparea_`jp'`shift' != 0
						lab var TC_in_act_jparea_`jp'`shift' "# of consecutive time periods with active AoE"	
					}
				local est_jp = subinstr(`"`estimation'"',"`evariable'","`complex'",1)
				local est_jp = subinstr(`"`est_jp'"',"`evariable'","in_act_jparea_`jp'`shift'",.)
				}

			if ("`complex'" == "") local est_jp = subinstr(`"`estimation'"',"`evariable'","in_act_jparea_`jp'`shift'",.)
			`est_jp'
			
			matrix b=e(b)
			matrix V=e(V)
			local indepvars : colnames b
			forvalues i = 1 / `n_b' {
				local colbase: word `i' of `indepvars' // `e(indepvars)'
				local colbase `=subinstr("`colbase'","_`jp'`shift'","",.)'
				if strpos("`colbase'","#") & strpos("`colbase'","c.") local col_int `=subinstr("`colbase'","c.","co.",.)'
				tokenize `aoe_bs_jp'
				forvalues j = 1 / `n_aoe_b' {	
					if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
						scalar b`j' = b[1,`i']
						scalar se`j' = sqrt(V[`i',`i'])
						if "`go_t'" == "1" scalar pv`j' = 2 * ttail(e(df_r),abs(b`j'/se`j'))
						if "`go_z'" == "1" scalar pv`j' = 2 * (1-normal(abs(b`j'/se`j')))
						}		
					}
				if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
					tokenize `x_bs'
					forvalues j = 1 / `n_x_b' {	
						forvalues k = 1 / `n_ovrep' { 
							if "`: word `j' of `x_bs''" == "`: word `k' of `ovreport''" {
								if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
									local x_lab_`j'_ev "`colbase'"
									scalar xb`j' = b[1,`i']
									di "Beta(`x_lab_`j'_ev') = "xb`j'
									scalar xse`j' = sqrt(V[`i',`i'])
									if "`go_t'" == "1" scalar xpv`j' = 2 * ttail(e(df_r),abs(xb`j'/xse`j'))
									if "`go_z'" == "1" scalar xpv`j' = 2 * (1-normal(abs(xb`j'/xse`j'))) //prob with logit & other var
									}					
								}
							}
						}	
					}	
				}	
			scalar p_jitter = `jp'
			foreach j of local diag {
				if !inlist("`j'","p_F","p_chi2") scalar `j' = e(`j')
				if "`j'" == "p_F" scalar `j' = Ftail(e(df_m),e(df_r),e(F))
				if "`j'" == "p_chi2" scalar `j' = chi2tail(e(df_m),e(chi2))
				}
			local cmd = e(cmd)
			post `memhold' (p_jitter) `tsb' `b_seqb' `se_seqb' `pv_seqb' `xb_seqb' `xse_seqb' `xpv_seqb' ///
				`diagb' ("`cmd'") ("`method'") ("`egridmatch'")
				
			if !inlist("`progress'","no","detail") nois _dots `s_dot' 0	
			if ("`progress'" == "detail") {
				local dot_return 0
				forvalues j = 1 / `n_aoe_b' {	
					if b`j' == 0  local dot_return 1
					}
				nois _dots `s_dot' `dot_return'
				}
			} 
		restore	
		}	
	postclose `memhold'	
	if ("`progress'" != "no") nois di " "			
	
********************************************************************************
	*** 1.3 PLACEBO AOE VISUALISATION ***

	use "`pos_pareas'", clear
	
		**************************************
		* Coordinates for added outline graph 
		gen x_ol = _X-x_cent if p_jitter == 0 
			lab var x_ol "longitude coordinates for outline graph"
		gen y_ol = _Y-y_cent if p_jitter == 0 
			lab var y_ol "latitude coordinates for outline graph"
			
		sum x_ol if p_jitter == 0
		local n_ol = r(N)	
		
		global xol ""
		global yol ""
		global olid ""
		
		levelsof aux if p_jitter == 0 & _X !=., local(olsteps) 
		foreach c in x y {
			foreach i of local olsteps {
				sum `c'_ol if p_jitter == 0 & aux == `i'
				global `c'ol $`c'ol `r(mean)'
				if "`c'" == "x" sum _ID if p_jitter == 0 & aux == `i'
				if "`c'" == "x" global olid $olid `r(mean)'
				}
			}
		if "`fixpoint'" == "anchored" {
			sum _ID if anchor == 1, meanonly 
			local anchor_id = `r(mean)'
			}
		**************************************
		
	keep if aux == 1
	tempfile jitter_add
	save `jitter_add'

	use "`pos_pout'", clear	
	
	merge m:1 p_jitter using "`jitter_add'", keepusing(x_adj y_adj)
		drop _merge

	* Note: If b = 0 & se = 0 the area is too small to include any incident, and thus, 
	*			the corresponding placebo variable was ommitted from the estimation.
	if ("`evreport'" == "" | "`complex'" == "") {
		local b b1
		local pv pv1
		local b_title `evariable'
		}
		
	forvalues j = 1 / `n_aoe_b' {	
		gen aux = 1 if b`j' == 0 & se`j' == 0 & missing(pv`j')
		qui replace b`j' = . if aux == 1
		qui replace se`j' = . if aux == 1
		lab var b`j' "b [`aoe_lab_`j'_ev']"
		lab var se`j' "se [`aoe_lab_`j'_ev']"
		lab var pv`j' "p-value [`aoe_lab_`j'_ev']"
		drop aux	
		if ("`evreport'" != "" & "`complex'" != "") {
			if "`evreport'" == "`aoe_lab_`j'_ev'" {
				local b b`j'
				local pv pv`j'
				local b_title `aoe_lab_`j'_ev'
				}
			}
		}
	if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
		forvalues j = 1 / `n_x_b' {	
			cap lab var xb`j' "b [`x_lab_`j'_ev']"
			cap lab var xse`j' "se [`x_lab_`j'_ev']"
			cap lab var xpv`j' "p-value [`x_lab_`j'_ev']"
			cap format xpv`j' %9.4f
			}
		}
	
	* Proportionality weight
	gen double `b'_abs = round(abs(`b'), 0.000001)	
		lab var `b'_abs "absolute proportionality weight"

	// Problem: no axes with weights when no observations
	forvalues scg = 1/6 {		
		local usew`scg' [aw = `b'_abs]
		}
	count if `b' < 0 & `pv' < 0.01		
		if r(N) == 0 local usew1 ""
	count if `b' < 0 & `pv' < 0.05 & `pv' >= 0.01 & !missing(`pv')
		if r(N) == 0 local usew2 ""
	count if `b' < 0 & `pv' < 0.1 & `pv' >= 0.05 & !missing(`pv')
		if r(N) == 0 local usew3 ""
	count if `b' > 0 & !missing(`pv') & `pv' < 0.01 
		if r(N) == 0 local usew4 ""
	count if `b' > 0 & !missing(`pv') & `pv' < 0.05 & `pv' >= 0.01 & !missing(`pv') 
		if r(N) == 0 local usew5 ""
	count if `b' > 0 & !missing(`pv') & `pv' < 0.1 & `pv' >= 0.05 & !missing(`pv')
		if r(N) == 0 local usew6 ""
	
	if "`templacebo'" != "" sort p_jitter ts
	
	// generate outline of real area
	d `b'
	if `r(N)' < `n_ol' { 
		set obs `n_ol'		// expand to match number of coordinates
		}
	gen ID_ol = "" //tmp
	foreach c in x y {
		tokenize $`c'ol
		gen `c'_ol = ""
		local `c'ol_n = wordcount("$`c'ol")
		forvalues ol = 1/``c'ol_n' {
			replace `c'_ol = "``ol''" in `ol'
			if "`c'" == "x" replace ID_ol = "`: word `ol' of $olid'" in `ol' //tmp
			}
			destring `c'_ol, replace
			if "`method'" == "geodetic" replace `c'_ol = `c'_ol * 100 // km adjustment
		}
	
	destring ID_ol, replace	
	local g_suba ""
	if "`multiarea'" == "" {
		local gl_suba (connected y_ol x_ol, mcolor(none) msymbol(none) lcolor(black`c_damp80') lwidth(medthick))
		local go_suba `"9 "observed area""'
		}
	if "`fixpoint'" == "anchored" {
		gen anchor = 1 if ID_ol == `anchor_id'
		sort anchor
		levelsof ID_ol, local(ID_ol)
		foreach suba of local ID_ol {
			local oll "lpattern(shortdash)"
			if `suba' == `anchor_id' local oll "lpattern(solid)"
			local gl_suba `gl_suba' (connected y_ol x_ol if ID_ol == `suba' , mcolor(none) msymbol(none) lcolor(black`c_damp80') `oll' lwidth(medthick))
			}
		local go_suba `"9 "observed (anchor) area" 10 "time-variant area(s)""'
		local gh_suba hole(10 11)
		}
	if inlist("`fixpoint'","joint","split") {
		levelsof ID_ol, local(ID_ol)
		foreach suba of local ID_ol {
			local gl_suba `gl_suba' (connected y_ol x_ol if ID_ol == `suba' , mcolor(none) msymbol(none) lcolor(black`c_damp80') lwidth(medthick))
			}
		local go_suba `"9 "time-variant areas""'
		local gh_suba hole(10 11)
		}		
		
	sort p_jitter
	save "`outres'/pos_pout`outname'.dta", replace	
	
	if "`templacebo'" == "" {
		if "`method'" == "geodetic" local jilab "in km" 
		if "`method'" == "degree" local jilab "in degree"
		nois graph twoway ///
			(scatter y_adj x_adj if `b' > 0 & !missing(`b') & `pv' < 0.01 `usew4', mcolor(`c_p1') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' > 0 & !missing(`b') & `pv' < 0.05 & `pv' >= 0.01 & !missing(`b') `usew5', mcolor(`c_p5') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' > 0 & !missing(`b') & `pv' < 0.1 & `pv' >= 0.05 & !missing(`b') `usew6', mcolor(`c_p10') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' < 0 & `pv' < 0.01 `usew1', mcolor(`c_n1') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' < 0 & `pv' < 0.05 & `pv' >= 0.01 & !missing(`b') `usew2', mcolor(`c_n5') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' < 0 & `pv' < 0.1 & `pv' >= 0.05 & !missing(`b') `usew3', mcolor(`c_n10') mlwidth(none)) ///
			(scatter y_adj x_adj if `b' == 0 | (`pv' >= 0.1 & !missing(`b')), mcolor(gs12) mlwidth(none)) ///	
			(scatter y_adj x_adj if missing(`pv'), m("X") mcolor(gs12)) `gl_suba' , ///	
			xtick(-`pos_r' (`pos_s') `pos_r') ytick(-`pos_r' (`pos_s') `pos_r') ///
			xtitle("longitude jitter (`jilab')") ytitle("latitude jitter (`jilab')") ///	
			xsize(4) ysize(4) ylabel(, nogrid) ///
			graphregion(color(white)) bgcolor(white) legend(order(1 "{&beta}{super:***}> 0" ///
			4 "{&beta}{super:***}< 0" 7 "{&beta}  =  0" 2 "{&beta}{super:**} > 0" ///
			5 "{&beta}{super:**} < 0" 8 "no incident variation" 3 "{&beta}{super:*}  > 0"  ///
			6 "{&beta}{super:*}  < 0" `go_suba') cols(3) `gh_suba') ///
			title("AoE position placebo estimates for {bf:`b_title'}", size(medium) color(black)) scale(0.75) ///
			name(pos_placebo, replace)
		}
		
	if "`templacebo'" != "" {
		local gr_stack ""		
		local gr_stack_in ""		
		if "`method'" == "geodetic" local jilab "Position shifts" "in km" 
		if "`method'" == "degree" local jilab "Position shifts" "in degree"
		
		gen gmat_col_pos = . 
		gen gmat_row_pos = . 

		// Position in matrix graph
		levelsof y_adj, local(col_shift)
		tokenize `col_shift'
		forvalues i = 1/ `=2*`pos_steps'+1' {
			replace gmat_col_pos = `i' + 1 if x_adj == ``i''
			replace gmat_row_pos = `i' if y_adj == ``i''*(-1)
			}
			
		local gmatpos_max = (2*`pos_steps'+1+1)^2

		// running number
		gen gmat_pos = (gmat_row_pos - 1) * (2 * `pos_steps' + 1 + 1) + gmat_col_pos

		// fake graphs for overview panel "labels": X
		levelsof x_adj, local(col_shift) //tmp
		tokenize `col_shift'
		forvalues i = 1 / `=2*`pos_steps'+1' {
			graph twoway (line `b' ts if cmd == "wtf" & ts == 0), ///
				text(0 0 "{&Delta}X = `=``i''`geomult''", size(*2.5) place(s) height(10)) yline(0, lcolor(black) lwidth(thick)) ///
				ylab(, nogrid nolabels noticks) ///
				ysize(2) ytitle("") yscale(lwidth(none)) xtitle("") xscale(lwidth(none)) ///
				graphregion(color(gs14)) plotr(color(gs14)) name(pos_placebo_ts_`=(2*`pos_steps'+1+1)*(2*`pos_steps'+1)+1+`i'', replace) nodraw
			local gr_stackx `gr_stackx' pos_placebo_ts_`=(2*`pos_steps'+1+1)*(2*`pos_steps'+1)+1+`i''
			}
			
		// fake graphs for overview panel "labels": Y
		levelsof y_adj, local(col_shift) //tmp
		tokenize `col_shift'
		forvalues i = 1 / `=2*`pos_steps'+1' {
			graph twoway (line `b' ts if cmd == "wtf" & ts == 0), ///
				text(0 0 "{&Delta}Y = `=``i''`geomult'*(-1)'", size(*2.5) place(w) width(20)) xline(0, lcolor(black) lwidth(thick)) ///
				ylab(, nogrid nolabels noticks) ///
				xsize(2) ytitle("") yscale(lwidth(none)) xtitle("") xscale(lwidth(none)) ///
				graphregion(color(gs14)) plotr(color(gs14)) name(pos_placebo_ts_`=1+(`i'-1)*(2*`pos_steps'+1+1)', replace) nodraw
			local gr_stacky `gr_stacky' pos_placebo_ts_`=1+(`i'-1)*(2*`pos_steps'+1+1)'
			}
		// fake graph for unit display	
		graph twoway (line `b' ts if cmd == "wtf" & ts == 0), ///
			text(0 0 "`jilab'", size(*2.5) place(w) width(20)) ///
			ylab(, nogrid nolabels noticks) ///
			ytitle("") yscale(lwidth(none)) xtitle("") xscale(lwidth(none)) ///
			graphregion(color(gs14)) plotr(color(gs14)) name(pos_placebo_ts_`=(2*`pos_steps'+1+1)*(2*`pos_steps'+1)+1', replace) nodraw	

		local gr_stack_in pos_placebo_ts_`=(2*`pos_steps'+1+1)*(2*`pos_steps'+1)+1'
		if ("`progress'" != "no") {
			nois di " "	
			nois di " "	
			nois di as text "Creating temporal AoE position placebo diagnostic graphs"
			}
		local gmat_in ""
		levelsof gmat_pos, local(gmat_in)
		foreach gmat of local gmat_in {
			sum x_adj if gmat_pos == `gmat', meanonly
			local dx = r(mean)
			sum y_adj if gmat_pos == `gmat', meanonly
			local dy = r(mean)	
			local plc white
			if `dy' == 0 & `dx' == 0 local plc gs14
			local xyline yline(0, lcolor(gs14)) xline(0, lcolor(black))
			local empty_X ""
			sum `b' if gmat_pos == `gmat', meanonly
			if r(N) == 0 {
				local xyline "" 
				local empty_X title("X", size(*3) color(gs14) ring(0) position(0))
				}
			graph twoway (bar `b' ts if gmat_pos == `gmat' & `b' > 0 & `pv' < 0.01 & inrange(ts,`sr_start',`sr_end'), color(`c_p1') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & `b' > 0 & `pv' < 0.05 & `pv' >= 0.01 & inrange(ts,`sr_start',`sr_end'), color(`c_p5') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & `b' > 0 & `pv' < 0.1 & `pv' >= 0.05 & inrange(ts,`sr_start',`sr_end'), color(`c_p10') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & `b' < 0 & `pv' < 0.01 & inrange(ts,`sr_start',`sr_end'), color(`c_n1') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & `b' < 0 & `pv' < 0.05 & `pv' >= 0.01 & inrange(ts,`sr_start',`sr_end'), color(`c_n5') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & `b' < 0 & `pv' < 0.1 & `pv' >= 0.05 & inrange(ts,`sr_start',`sr_end'), color(`c_n10') lcolor(white)) ///
				(bar `b' ts if gmat_pos == `gmat' & (`b' == 0 | (`pv' >= 0.05 & !missing(`pv'))) & inrange(ts,`sr_start',`sr_end'), color(gs12) lcolor(white)), ///
				`xyline' `empty_X' ylab(, nogrid nolabels noticks) yscale(lwidth(none)) xlab("") xscale(lwidth(none)) /// 
				ytitle("") xtitle("") legend(off) ///
				graphregion(color(gs14)) plotr(color(`plc')) bgcolor(white) name(pos_placebo_ts_`gmat', replace) nodraw
			local gr_stack_in `gr_stack_in' pos_placebo_ts_`gmat'
			}

		local gr_stack `gr_stacky' `gr_stackx' `gr_stack_in'

		local gr_stack2 ""
		local gmat_in ""
		forvalues j = 1 / `gmatpos_max' {
			local check : word `j' of `gr_stack'
			foreach x of local gr_stack {
				local run `=subinstr("`x'","pos_placebo_ts_","",1)'
				if "`j'" == "`run'" local gr_stack2 `gr_stack2' pos_placebo_ts_`run'
				if "`j'" == "`run'" local gmat_in `gmat_in' `run'
				}
			}
			
		local gmat_holes ""
		local gmat_list ""
		forvalues j = 1 / `gmatpos_max' {
			local gmat_list `gmat_list' `j'
			}

		local gmat_holes: list gmat_list - gmat_in
		if ("`progress'" != "no") nois di as text "Combining temporal AoE position placebo diagnostics graphs"
		graph combine `gr_stack2', graphregion(color(gs14) margin(l=22 r=22)) ///
			cols(`=2*`pos_steps'+1+1') imargin(0 0 0 0) holes(`gmat_holes') ycommon ///
			title("AoE position placebo estimates for {bf:`b_title'}", size(medium) color(black))
		}		
		
		graph save "`outvis'/pos_pdiag`outname'.gph", replace	
		graph export "`outvis'/pos_pdiag`outname'.png", replace
		}
	

********************************************************************************
*** 2 ROTATION PLACEBO TEST													 ***	
********************************************************************************

if ("`rotation'" != "") {
	if ("`method'" == "degree") {
		
		if ("`progress'" != "no") {
			nois di " "
			*if ("`progress'" != "no") nois di in red "*** ROTATION AOE PLACEBOS ***"	
			nois di "{opt Placebo type}:      {opt Rotation AOE placebos}"
			nois di " "
			}
			
		*** 2.1 PLACEBO AREA GENERATION ***
		
		clear
		
		use "`areadata'", clear
		
			* calculate the centroid
			gen aux=_n
			
			if inlist("`fixpoint'","anchored","split") local centgen "bysort _ID  `ea_start' `ea_end' (aux):"

			* calculate the centroid
			`centgen' gen a=0.5*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen A=total(a)
				replace A=abs(A)
			`centgen' gen c_x=(_X+_X[_n+1])*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen x_cent=total(c_x)
				replace x_cent=abs(1/(6*A)*x_cent)
			`centgen' gen c_y=(_Y+_Y[_n+1])*(_Y*_X[_n+1]-_Y[_n+1]*_X)
			`centgen' egen y_cent=total(c_y)
			replace y_cent=abs(1/(6*A)*y_cent)	
			
		if "`fixpoint'" == "anchored" {
			* Set anchor centroid (earliest area)
			sum `ea_start'
			gen anchor = 1 if `ea_start' == `r(min)'
			sum x_cent if anchor == 1, meanonly 
			replace x_cent = `r(mean)'
			sum y_cent if anchor == 1, meanonly 
			replace y_cent = `r(mean)'
			local kp anchor
			}
			
		keep _ID _Y _X ?_cent aux `kp'
		
		* Derive recentred vertices (centroid is assumed to be 0/0))
		gen x_c = _X-x_cent
		gen y_c = _Y-y_cent

		* Rotate for each of the rotation steps
		local rot_steps = floor((`rot_b')/`rot_s')	//@rotmod
		forvalues rot = 1/`rot_steps' {
			gen rotation`rot'_pos = `rot'*`rot_s'
			gen rotation`rot'_neg = -`rot'*`rot_s'
			
			gen _X`rot'_pos = -y_c*sin(-`rot'*`rot_s'*_pi/180)+x_c*cos(-`rot'*`rot_s'*_pi/180)
				replace _X`rot'_pos = _X`rot'_pos+x_cent
			gen _X`rot'_neg = -y_c*sin(`rot'*`rot_s'*_pi/180)+x_c*cos(`rot'*`rot_s'*_pi/180)
				replace _X`rot'_neg = _X`rot'_neg+x_cent				
				
			gen _Y`rot'_pos = y_c*cos(-`rot'*`rot_s'*_pi/180)+x_c*sin(-`rot'*`rot_s'*_pi/180)
				replace _Y`rot'_pos = _Y`rot'_pos+y_cent
			gen _Y`rot'_neg = y_c*cos(`rot'*`rot_s'*_pi/180)+x_c*sin(`rot'*`rot_s'*_pi/180)
				replace _Y`rot'_neg = _Y`rot'_neg+y_cent				
			}

		rename _X _X0
		rename _Y _Y0
		gen rotation0 = 0
			
		reshape long _X _Y rotation, i(_ID aux) string
		
		drop ?_c
		egen p_rot = group(_j)
		drop _j

		sort p_rot _ID aux	//@mod
		
		gen orig=1 if rotation==0
			replace orig=0 if orig==.
			
		tempfile rot_pareas
		save `rot_pareas'
			save "`outpoly'/rot_pareas`outname'.dta", replace
			
		* Export to temporary placebo area reference files
		sum p_rot
		local rp_min = r(min)
		local rp_max = r(max)
		
		if ("`progress'" != "no") nois _dots 0, title(Generating rotation placebo areas (Total: `rp_max'))		
		forvalues rp = `rp_min'/`rp_max' {
			if ("`progress'" != "no") nois _dots `=`rp'+1' 0
			sum rotation if p_rot == `rp'
			local rot_`rp' = r(mean)
			tempfile rot_parea_`rp'
			savesome _all using "`rot_parea_`rp''" if p_rot == `rp'
			}
		}	

	if ("`method'" == "geodetic") {
		use "`area_geo'", clear

		if ("`progress'" != "no") {
			nois di " "	
			*if ("`progress'" != "no") nois di in red "*** ROTATION AOE PLACEBOS: {opt Geodetic} adjustments ***"	
			if ("`progress'" != "no") nois di "{opt Placebo type}:      {opt Rotation AOE placebos}"
			nois di " "	
			}
			
		local rot_steps = floor(`rot_b'/`rot_s')	//@rotmod
		
		if ("`progress'" != "no") nois _dots 0, title(Generating rotation placebo areas (Total: `=`rot_steps'*2+1'))	
		local count_pos = 0
		
		// Add random rotation component to original radian
		gen double rot_vertex0 = 0*`rot_s' + rad_vertex
		if ("`progress'" != "no") {
			local count_pos = `count_pos' + 1
			nois _dots `count_pos' 0	
			}
			
		forvalues rot = 1/`rot_steps' {
			if ("`progress'" != "no") {
				forvalues ct = 1 / 2 {
					local count_pos = `count_pos' + 1
					nois _dots `count_pos' 0 
					}
				}			
			gen double rot_vertex`rot'_pos = `rot'*`rot_s'*c(pi)/180 + rad_vertex
			gen double rot_vertex`rot'_neg = -`rot'*`rot_s'*c(pi)/180 + rad_vertex
			}
			
		reshape long rot_vertex, i(aux) string	
		split _j, p("_")
		destring _j1, replace
		replace _j1 = _j1*(-1) if _j2=="neg"
		rename _j1 rot_step
		drop _j _j2
		
		sort rot_step aux
		
		*** 2.2 ROTATING: Get new vertices in relation to jittered centroid & random rotation ***
		tempname d2r a b f
		scalar `d2r' = c(pi) / 180
		scalar `a' = 6378137
		scalar `f' = 1/298.257223563 // WGS-84 ellipsiod	
		scalar `b' = `a' - `a' * `f'
	
		local maxiter 1000
		
		local lon1 x_cent
		local lat1 y_cent

		// at least one lat/lon is a variable; others can be scalar or number
		tempvar touse
		gen `touse' = 1 `if' `in'
		qui replace `touse' = 0 if mi(`lat1',`lon1')
		qui replace `touse' = 0 if `touse' == .					
		
		tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
			cos_Sqalpha uSq A B 
			
		gen double `sin_alpha1'	= sin(rot_vertex)
		gen double `cos_alpha1' = cos(rot_vertex)
		gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
		gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
		gen double `sin_U1' = `tan_U1'*`cos_U1' 
		gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
		gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
		gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
		gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
		gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
		gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
		
		// iterations
		tempvar cont
		gen `cont' = `touse'	
		foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
			tempvar `v'
			gen double ``v'' = .
			}
		qui replace `sigma' = disto*1000/(`b'*`A')									
		
		local iter 0
		local more 1
		while `++iter' < `maxiter' & `more' {
			qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
			qui replace `sin_sigma' = sin(`sigma') if `cont'					
			qui replace `cos_sigma' = cos(`sigma') if `cont'					
			qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
				`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
			qui replace `sigma_old' = `sigma' if `cont'
			qui replace `sigma' = disto*1000/(`b'*`A') + `delta_sigma' if `cont'	

			// mark out observations that have converged
			qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
			// we are done if all observations have converged
			sum `cont', meanonly
			local more = r(max)
			}		
		
		tempvar x lambda C L lambda2

		gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
		gen double _Y_r_geo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
		gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
		gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
		gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
		gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
		gen double _X_r_geo = `lambda2'*180/c(pi)
		qui replace _Y_r_geo = _Y_r_geo*180/c(pi)
		
		drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
			`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
			`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont' `touse'
		cap drop __0000*

		bysort rot_step: gen auxj = _n 
		bysort auxj (rot_step): gen auxj2 = _n if auxj == 1
		bysort rot_step: egen p_rot = max(auxj2)
		
		gen orig=1 if rot_step == 0 
			replace orig=0 if orig==.	
		
		sort p_rot aux	
		drop auxj auxj2 _Y _X
		rename _Y_r_geo _Y
		rename _X_r_geo _X
		
		sum p_rot
		local rp_min = r(min)
		local rp_max = r(max)
		
		gen rotation = rot_step * `rot_s'

		tempfile rot_pareas
		save `rot_pareas'
			save "`outpoly'/rot_pareas`outname'.dta", replace // tmp	

		forvalues rp = `rp_min'/`rp_max' {
			tempfile rot_parea_`rp'
			savesome _all using	`rot_parea_`rp'' if p_rot == `rp'
			}		
		} // END geo
		
********************************************************************************		
	
	* Identify incidents in various areas
	if ("`progress'" != "no") {
		nois di " "					
		nois di " "	
		nois _dots 0, title(Identifying cells in rotation placebo areas via `egridmatch' link (Total: `=`rp_max''))		
		}
		
	if "`egridmatch'" == "centroid" {	
		use "`celldata'", clear
		
		forvalues rp = `rp_min'/`rp_max' {
			if ("`progress'" != "no") nois _dots `=`rp'+1' 0
			if "`multiarea'" == "" {
				geoinpoly y_cent x_cent using "`rot_parea_`rp''", unique 	
				rename _ID in_rparea_`rp'										
				qui replace in_rparea_`rp'=1 if !missing(in_rparea_`rp')
				qui replace in_rparea_`rp'=0 if missing(in_rparea_`rp')
				gen in_act_rparea_`rp'=(in_rparea_`rp'==1 & `etime'==1) if !missing(`etime')
				}
			if "`multiarea'" != "" {
				geoinpoly y_cent x_cent using "`rot_parea_`rp''"	
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
				replace _ID = 1 if inrange(`tid',ea_start,ea_end) & !missing(_ID)
				replace _ID = 0 if !inrange(`tid',ea_start,ea_end) | missing(_ID)
				bysort `grid' `tid': egen in_act_rparea_`rp' = max(_ID)
				drop ea_* _ID
				duplicates drop `grid' `tid' in_act_rparea_`rp', force
				}
			lab var in_act_rparea_`rp' "Centroid of (grid) cell in active operation area"
			}
		}
			
	if "`egridmatch'" == "vertex" {		
		if "`multiarea'" == "" {
			use "`cellpoly'", clear
		forvalues rp = `rp_min'/`rp_max' {
				if ("`progress'" != "no") nois _dots `=`rp'+1' 0
						
					geoinpoly _Y _X using "`rot_parea_`rp''", unique
					rename _ID aux	
					qui replace aux=1 if !missing(aux)
					qui replace aux=0 if missing(aux)
					bysort `grid': egen in_rparea_`rp'=max(aux)
					drop aux
					}
				
			drop _X _Y
			duplicates drop

			merge 1:m `grid' using "`celldata'"
				keep if _merge==3
				drop _merge
			
			forvalues rp = `rp_min'/`rp_max' {
				gen in_act_rparea_`rp'=(in_rparea_`rp'==1 & `etime'==1) if !missing(`etime')
				lab var in_act_rparea_`rp' "Vertex of (grid) cell in active operation area"
				}
			}	
		
		if "`multiarea'" != "" {
			forvalues rp = `rp_min'/`rp_max' {
				use "`cellpoly'", clear
				if ("`progress'" != "no") nois _dots `=`rp'+1' 0
				geoinpoly _Y _X using "`rot_parea_`rp''"
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
	
				rename _ID aux	
				qui replace aux=0 if missing(aux)
		
				drop if missing(_X)
				drop _X _Y
				duplicates drop
		
				reshape wide ea_start ea_end, i(`grid') j(aux)
				
				merge 1:m `grid' using "`celldata'", keepusing(`tid')
	
				gen in_act_rparea_`rp' = 0
					lab var in_act_rparea_`rp' "Vertex of (grid) cell in active operation area"

				foreach sa of local subareas {
					replace	in_act_rparea_`rp' = 1 if inrange(`tid',ea_start`sa',ea_end`sa') ///
						& !missing(ea_start`sa') & !missing(ea_end`sa')	
					}	
				drop ea_start* ea_end* _merge
				tempfile mat_`rp'
				save `mat_`rp''
				}	
			use "`celldata'", clear
			forvalues rp = `rp_min'/`rp_max' {
				merge 1:1 `grid' `tid' using `mat_`rp''
					drop _merge
				}
			}
		}		
		
	tempfile rot_pest
	save `rot_pest'	
	save "`outres'/rot_pest`outname'.dta", replace
		
	* Optional / Test mapping
	* NOTE: polygon & data can't handle temp files (messed modification)
	if ("`mapdb'" != "") & ("`mapco'" != "") {
		use "`mapdb'", clear
		spmap using "`mapco'", ///
			id(_ID) fcolor(gs12) ocolor(gs12) ///
			polygon(data("`outpoly'/rot_pareas`outname'.dta") by(orig) osize(vthin) ocolor(black cranberry)) ///
			point(data("`outpoly'/rot_pareas`outname'.dta") x(x_cent) y(y_cent) ///
			size(vsmall) fcolor(cranberry))
		graph save "`outvis'/rot_pareas`outname'.gph", replace	
		graph export "`outvis'/rot_pareas`outname'.png", replace	
		}	
	
********************************************************************************	
	*** 2.2 PLACEBO AOE ESTIMATIONS ***
		
	tempname memhold
	tempfile rot_pout
	cap postclose `memhold'
	postfile `memhold' rotation `ts' `b_seq' `se_seq' `pv_seq' `xb_seq' `xse_seq' `xpv_seq' ///
		`diag' str30(cmd method egridmatch) using "`rot_pout'"
		
	use "`rot_pest'", clear
	
	if ("`progress'" != "no") nois di " "					
	if ("`progress'" != "no") nois di " "					
	if "`templacebo'" != "" {
		if ("`progress'" != "no") nois _dots 0, title(Estimation for rotation placebo areas ///
		(Total: `=`rp_max'' areas x `n_tprange' temp. shifts = `=(`rp_max')*`n_tprange''))
		}
	if "`templacebo'" == ""	{
		if ("`progress'" != "no") nois _dots 0, title(Estimation for rotation placebo areas (Total: `=`rp_max''))
		}
		
	local s_dot = 0
	forvalues rp = `rp_min'/`rp_max' {
		preserve
			use "`rot_pareas'", clear
			qui sum rotation if p_rot == `rp', meanonly
			local rot_`rp' = r(mean)
		restore	
		
		preserve
		
		if "`templacebo'" != "" {
			xtset `grid' `tid'

			tempvar aux_treated aux_t1 aux_tT t1 tT

			bysort `grid': egen `aux_treated' = max(in_act_rparea_`rp')

			sum `tid'
			local tmin = r(min)
			local tmax = r(max)
			gen `aux_t1' = in_act_rparea_`rp' if `tid' == `tmin'
			gen `aux_tT' = in_act_rparea_`rp' if `tid' == `tmax'
			bysort `grid': egen `t1' = max(`aux_t1')
			bysort `grid': egen `tT' = max(`aux_tT')
			
			local shift_suffix ""
			forvalues s = `s_start' / `s_end' {
				if `s' < 0 {
					gen in_act_rparea_`rp'_sn_`=`s' * (-1)' = F`=`s' * (-1)'.in_act_rparea_`rp'
					local shift_suffix `shift_suffix' _sn_`=`s' * (-1)'
					if "`seriesmargin'" == "zero" & !missing(in_act_rparea_`rp') ///
						replace in_act_rparea_`rp'_sn_`=`s' * (-1)' = 0 if missing(in_act_rparea_`rp'_sn_`=`s' * (-1)')
					if "`seriesmargin'" == "preserve" & !missing(in_act_rparea_`rp') {
						replace in_act_rparea_`rp'_sn_`=`s' * (-1)' = 1 if missing(in_act_rparea_`rp'_sn_`=`s' * (-1)') & `aux_treated' == 1 & `tT' == 1
						replace in_act_rparea_`rp'_sn_`=`s' * (-1)' = 0 if missing(in_act_rparea_`rp'_sn_`=`s' * (-1)') & (`aux_treated' == 0 | `tT' == 0)
						}
					}
				if `s' > 0 {
					gen in_act_rparea_`rp'_sp_`=`s'' = L`s'.in_act_rparea_`rp'
					local shift_suffix `shift_suffix' _sp_`=`s''
					if "`seriesmargin'" == "zero" & !missing(in_act_rparea_`rp') replace in_act_rparea_`rp'_sp_`=`s'' = 0 if missing(in_act_rparea_`rp'_sp_`=`s'')
					if "`seriesmargin'" == "preserve" & !missing(in_act_rparea_`rp') {
						replace in_act_rparea_`rp'_sp_`=`s'' = 1 if missing(in_act_rparea_`rp'_sp_`=`s'') & `aux_treated' == 1 & `t1' == 1
						replace in_act_rparea_`rp'_sp_`=`s'' = 0 if missing(in_act_rparea_`rp'_sp_`=`s'') & (`aux_treated' == 0 | `t1' == 0)
						}
					}
				}
			}		
		
		foreach shift in "" `shift_suffix' {
			local s_dot = `s_dot' + 1
			if "`shift'" == "" scalar ts = 0
			if strpos("`shift'","_sn_") scalar ts = `=subinstr("`shift'","_sn_","-",1)'
			if strpos("`shift'","_sp_") scalar ts = `=subinstr("`shift'","_sp_","",1)'

			if ("`complex'" != "") {
				xtset `grid' `tid'
				
				if strpos("`complex'","TS_`evariable'") {
					bysort `grid' (`tid'): gen TS_in_act_rparea_`rp'`shift' = sum(in_act_rparea_`rp'`shift') if in_act_rparea_`rp'`shift'==1
						replace TS_in_act_rparea_`rp'`shift' = 0 if TS_in_act_rparea_`rp'`shift'==. & !missing(in_act_rparea_`rp'`shift')		
						lab var TS_in_act_rparea_`rp'`shift' "# of cumulative time periods with active AoE"	
					}
				if strpos("`complex'","TC_`evariable'") {
					bysort `grid' (`tid'): gen TC_in_act_rparea_`rp'`shift' = sum(in_act_rparea_`rp'`shift') if in_act_rparea_`rp'`shift'==1
					replace TC_in_act_rparea_`rp'`shift' = 0 if in_act_rparea_`rp'`shift' == 0
					tempvar aux
					gen `aux' = TC_in_act_rparea_`rp'`shift'
						replace `aux' = L1.`aux' if `aux' == 0 & L1.`aux' != .
						replace `aux' = L1.`aux' if `aux' == TC_in_act_rparea_`rp'`shift' & L1.`aux' != .
						replace TC_in_act_rparea_`rp'`shift' = TC_in_act_rparea_`rp'`shift' - `aux' if TC_in_act_rparea_`rp'`shift' != 0
						lab var TC_in_act_rparea_`rp'`shift' "# of consecutive time periods with active AoE"	
					}
				local est_rp = subinstr(`"`estimation'"',"`evariable'","`complex'",1)
				local est_rp = subinstr(`"`est_rp'"',"`evariable'","in_act_rparea_`rp'`shift'",.)
				}

			if ("`complex'" == "") local est_rp = subinstr(`"`estimation'"',"`evariable'","in_act_rparea_`rp'`shift'",.)
			`est_rp'		

			matrix b=e(b)
			matrix V=e(V)
			local indepvars : colnames b
			forvalues i = 1 / `n_b' {
				local colbase: word `i' of `indepvars' //`e(indepvars)'
				local colbase `=subinstr("`colbase'","_`rp'`shift'","",.)'
				if strpos("`colbase'","#") & strpos("`colbase'","c.") local col_int `=subinstr("`colbase'","c.","co.",.)'
				tokenize `aoe_bs_rp'
				forvalues j = 1 / `n_aoe_b' {	
					if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
						scalar b`j' = b[1,`i']
						scalar se`j' = sqrt(V[`i',`i'])
						if "`go_t'" == "1" scalar pv`j' = 2 * ttail(e(df_r),abs(b`j'/se`j'))
						if "`go_z'" == "1" scalar pv`j' = 2 * (1-normal(abs(b`j'/se`j')))
						}		
					}
				if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
					tokenize `x_bs'
					forvalues j = 1 / `n_x_b' {	
						forvalues k = 1 / `n_ovrep' { 
							if "`: word `j' of `x_bs''" == "`: word `k' of `ovreport''" {
								if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
									local x_lab_`j'_ev "`colbase'"
									scalar xb`j' = b[1,`i']
									di "Beta(`x_lab_`j'_ev') = "xb`j'
									scalar xse`j' = sqrt(V[`i',`i'])
									if "`go_t'" == "1" scalar xpv`j' = 2 * ttail(e(df_r),abs(xb`j'/xse`j'))
									if "`go_z'" == "1" scalar xpv`j' = 2 * (1-normal(abs(xb`j'/xse`j'))) //prob with logit & other var
									}					
								}
							}
						}	
					}	
				}		
			scalar rotation = `rot_`rp''
			foreach j of local diag {
					if !inlist("`j'","p_F","p_chi2") scalar `j' = e(`j')
					if "`j'" == "p_F" scalar `j' = Ftail(e(df_m),e(df_r),e(F))
					if "`j'" == "p_chi2" scalar `j' = chi2tail(e(df_m),e(chi2))
					}	
			local cmd = e(cmd)
			post `memhold' (rotation) `tsb' `b_seqb' `se_seqb' `pv_seqb' `xb_seqb' `xse_seqb' `xpv_seqb' ///
				`diagb' ("`cmd'") ("`method'") ("`egridmatch'")
			
			if !inlist("`progress'","no","detail") nois _dots `s_dot' 0	
			if ("`progress'" == "detail") {
				local dot_return 0
				forvalues j = 1 / `n_aoe_b' {	
					if b`j' == 0  local dot_return 1
					}
				nois _dots `s_dot' `dot_return'
				}
			} 
		restore	
		}	
	postclose `memhold'
	if ("`progress'" != "no") nois di " "					

********************************************************************************	
	*** 2.3 PLACEBO AOE VISUALISATION ***
	
	use "`rot_pout'", clear
	
	// Note: If b = 0 & se = 0 the area is too small to include any incident, and thus, 
	//			the corresponding placebo variable was ommitted from the estimation.
	if ("`evreport'" == "" | "`complex'" == "") {
		local b b1
		local se se1
		local pv pv1
		local b_title `evariable'
		}
		
	forvalues j = 1 / `n_aoe_b' {	
		gen aux = 1 if b`j' == 0 & se`j' == 0 & missing(pv`j')
		qui replace b`j' = . if aux == 1
		qui replace se`j' = . if aux == 1
		lab var b`j' "b [`aoe_lab_`j'_ev']"
		lab var se`j' "se [`aoe_lab_`j'_ev']"
		lab var pv`j' "p-value [`aoe_lab_`j'_ev']"
		drop aux	
		if ("`evreport'" != "" & "`complex'" != "") {
			if "`evreport'" == "`aoe_lab_`j'_ev'" {
				local b b`j'
				local se se`j'
				local pv pv`j'
				local b_title `aoe_lab_`j'_ev'
				}
			}
		}
	if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
		forvalues j = 1 / `n_x_b' {	
			cap lab var xb`j' "b [`x_lab_`j'_ev']"
			cap lab var xse`j' "se [`x_lab_`j'_ev']"
			cap lab var xpv`j' "p-value [`x_lab_`j'_ev']"
			cap format xpv`j' %9.4f
			}
		}
		
	* Centre original zero degree angle
	replace rotation = rotation - 360 if rotation > 180 & !missing(rotation)
	if "`method'" == "degree" {
		expand 2 if rotation == 180, gen(exp)		//@rotmod
		replace rotation = -180 if exp == 1
		drop exp
		}
	
	gen `b'_cl95 = `b' - invttail(df_r,(1-95/100)/2) * `se'
	gen `b'_cu95 = `b' + invttail(df_r,(1-95/100)/2) * `se'

	sort rotation
	
		save "`outres'/rot_pout`outname'.dta", replace	
		
	if "`templacebo'" == "" {
		nois graph twoway (rarea `b'_cl95 `b'_cu95 rotation, color(gs14`c_damp80')) (line `b' rotation, color(black)),  ///
			yline(0) xline(0) ylabel(, nogrid) xlab(-`rot_b' (`rot_s') `rot_b', alt) /*xtick(-180 (`rot_s') 180)*/  /// 
			xtitle("rotation angle") legend(order(2 "{&beta}" 1 "95% CI")) ///
			title("AoE rotation placebo estimates for {bf:`b_title'}", size(medium) color(black)) ///
			graphregion(color(white)) bgcolor(white) name(rot_placebo, replace)
		}
	
	if "`templacebo'" != "" {
		local gr_stack ""
		tempvar y
		gen `y' = 0	
		levelsof ts if inrange(ts,`sr_start',`sr_end'), local(shift)
		if ("`progress'" != "no")  {
			nois di " "
			nois di " "
			nois di as text "Creating temporal AoE rotation placebo diagnostic graphs"
			}
		if "`tpgtype'" == "line" {
			foreach s of local shift {
				if `s' < 0 {
					local sadj n
					local s1 *(-1)
					}
				if `s' > 0 {
					local sadj p
					local s1 ""
					}
				local plc white
				local ylc gs14
				if `s' == 0 {
					local plc gs14
					local ylc white
					}
				graph twoway (line `y' rotation, color(cranberry) lwidth(vthin)) ///
					(rarea `b'_cl95 `b'_cu95 rotation if ts == `s', color(`ylc'`c_damp80')) ///
					(line `b' rotation if ts == `s', color(black)) ///
					(line `b' rotation if cmd == "wtf", yaxis(2)) , ///
					xline(0, lcolor(cranberry)) ///
					ylab(, nogrid /*nolabels noticks*/ alternate) /*yscale(lwidth(none))*/ ///
					ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal)) xlab("") xscale(lwidth(none)) /// 
					ytitle("") xtitle("") legend(off) ///
					graphregion(color(`plc')) plotr(color(`plc')) bgcolor(gs14) name(rot_placebo_ts_`sadj'`=`s'`s1'', replace) nodraw
				local gr_stack `gr_stack' rot_placebo_ts_`sadj'`=`s'`s1''
				}
			// fake rotation label graph
			graph twoway (line `y' rotation, lwidth(none)) ///
				(line `b' rotation if cmd == "wtf", yaxis(2)), ///
				ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal) color(gs14)) ///
				ylabel(, nogrid /*nolabel*/ tlcolor(gs14) labcolor(gs14) alternate) yscale(lwidth(none)) ytitle("") //
				xlab(-`rot_b' (`rot_s') `rot_b', alt) /*xtick(-180 (`rot_s') 180)*/  /// 
				yscale(axis(2) lwidth(none)) xtitle("rotation angle") legend(off) ///
				graphregion(color(gs14)) plotr(color(gs14)) bgcolor(white) name(rot_placebo_ts_lab, replace) nodraw
				}
			
		if "`tpgtype'" != "line" {
			foreach s of local shift {
				if `s' < 0 {
					local sadj n
					local s1 *(-1)
					}
				if `s' > 0 {
					local sadj p
					local s1 ""
					}
				local plc white
				local ylc gs14
				if `s' == 0 {
					local plc gs14
					local ylc white
					}
	
				graph twoway (bar `b' rotation if ts == `s' & `b' > 0 & `pv' < 0.01, color(`c_p1') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & `b' > 0 & `pv' < 0.05 & `pv' >= 0.01, color(`c_p5') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & `b' > 0 & `pv' < 0.1 & `pv' >= 0.05, color(`c_p10') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & `b' < 0 & `pv' < 0.01, color(`c_n1') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & `b' < 0 & `pv' < 0.05, color(`c_n5') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & `b' < 0 & `pv' < 0.1 & `pv' >= 0.05, color(`c_n10') lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(bar `b' rotation if ts == `s' & (`b' == 0 | (`pv' >= 0.05 & !missing(`pv'))), color(gs12) lcolor(white) barwidth(`=0.9*`rot_s'') lwidth(none)) ///
					(line `b' rotation if cmd == "wtf", yaxis(2))  ///
					(line `y' rotation, color(`ylc') lwidth(thin)), ///
					xline(0, lcolor(black)) ///
					ylab(, nogrid /*nolabels noticks*/ alternate) ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal)) ///
					xlab("") xscale(range(-`rot_b' `rot_b') lwidth(none)) ytitle("") xtitle("") legend(off) ///
					graphregion(color(`plc')) plotr(color(`plc')) bgcolor(gs14) name(rot_placebo_ts_`sadj'`=`s'`s1'', replace) nodraw
				local gr_stack `gr_stack' rot_placebo_ts_`sadj'`=`s'`s1''
				}
			// fake rotation label graph
			graph twoway (bar `y' rotation, barwidth(`=0.9*`rot_s'') lwidth(none)) ///
				(line `b' rotation if cmd == "wtf", yaxis(2)), ///
				ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal) color(gs14)) ///
				ylabel(, nogrid tlcolor(gs14) labcolor(gs14) alternate) yscale(lwidth(none)) ytitle("") ///
				xlab(-`rot_b' (`rot_s') `rot_b', alt) /*xtick(-180 (`rot_s') 180)*/  /// 
				yscale(axis(2) lwidth(none)) xtitle("rotation angle") legend(off) ///
				graphregion(color(gs14)) plotr(color(gs14)) bgcolor(white) name(rot_placebo_ts_lab, replace) nodraw
			}
		local n_gr: word count `gr_stack'
		local gr_stack_in ""		
		tokenize `gr_stack'
		forvalues i = 1 / `n_gr' {
			local gr_stack_in `gr_stack_in' ``=`n_gr'-`i'+1''
			}
		if ("`progress'" != "no") nois di as text "Combining temporal AoE rotation placebo diagnostics graphs"
		graph combine `gr_stack_in' rot_placebo_ts_lab, graphregion(color(gs14)) ///
			title("AoE rotation placebo estimates for {bf:`b_title'}", size(medium) color(black)) ///	
			cols(1) imargin(0 0 0 0) ycommon		
		}	
		
		graph save "`outvis'/rot_pdiag`outname'.gph", replace	
		graph export "`outvis'/rot_pdiag`outname'.png", replace
	}
	}
	
********************************************************************************
*** 3 SCALE PLACEBO TEST													 ***	
********************************************************************************

qui {
if ("`scale'" != "") {

	if ("`progress'" != "no") {
		nois di " "	
		*if ("`progress'" != "no") nois di in red "*** SCALE AOE PLACEBOS ***"	
		if ("`progress'" != "no") nois di "{opt Placebo type}:      {opt Scale AOE placebos}"
		nois di " "
		}
		
	*** 3.1 PLACEBO AREA GENERATION ***
	
	if ("`method'" == "degree") {

		clear
		
		use "`areadata'", clear
		
			* calculate the centroid
			gen aux=_n
			
			if inlist("`fixpoint'","anchored","split") local centgen "bysort _ID  `ea_start' `ea_end' (aux):"

			* calculate the centroid
			`centgen' gen a=0.5*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen A=total(a)
				replace A=abs(A)
			`centgen' gen c_x=(_X+_X[_n+1])*(_X*_Y[_n+1]-_X[_n+1]*_Y)
			`centgen' egen x_cent=total(c_x)
				replace x_cent=abs(1/(6*A)*x_cent)
			`centgen' gen c_y=(_Y+_Y[_n+1])*(_Y*_X[_n+1]-_Y[_n+1]*_X)
			`centgen' egen y_cent=total(c_y)
			replace y_cent=abs(1/(6*A)*y_cent)	
			
		if "`fixpoint'" == "anchored" {
			* Set anchor centroid (earliest area)
			sum `ea_start'
			gen anchor = 1 if `ea_start' == `r(min)'
			sum x_cent if anchor == 1, meanonly 
			replace x_cent = `r(mean)'
			sum y_cent if anchor == 1, meanonly 
			replace y_cent = `r(mean)'
			local kp anchor
			}
			
		keep _ID _Y _X ?_cent aux `kp'
		
			* Derive vectors from centroid to vertices
			gen x_vec = _X-x_cent
			gen y_vec = _Y-y_cent
				
			* Retrieve closest starting point to zero
			if "`scale_lb'" == "" local scale_lb = 1-(1/`scale_s'-1) * `scale_s'
			local s = 0
			foreach elast of numlist `scale_lb' (`scale_s') `scale_ub' {
				local s = `s' + 1
				gen scale`s' = `elast'
				gen _X`s' = (`elast') * x_vec + x_cent
				gen _Y`s' = (`elast') * y_vec + y_cent
				}
				
			drop _X _Y

			reshape long _X _Y scale, i(_ID aux)

			drop /*_ID*/ *_vec
			*gen _ID = _j
			rename _j p_scale
			qui replace p_scale = 0 if scale == 1
			qui replace p_scale = p_scale - 1 if scale > 1 & !missing(scale)
		
		sort p_scale _ID aux
		}
		
	if ("`method'" == "geodetic") {
		use "`area_geo'", clear

		local scale_lb = 1-(1/`scale_s'-1) * `scale_s'
		local s = 0
		foreach elast of numlist `scale_lb' (`scale_s') `scale_ub' {
			local s = `s' + 1
			gen scale`s' = `elast'
			}
		
		reshape long scale, i(aux)
		
		// sort scale aux
		tempvar dist_shift
		gen double `dist_shift' = 1000 * disto * scale

		tempname d2r a b f
		scalar `d2r' = c(pi) / 180
		scalar `a' = 6378137
		scalar `f' = 1/298.257223563 // WGS-84 ellipsiod	
		scalar `b' = `a' - `a' * `f'
		
		local maxiter 1000
	
		local lon1 x_cent
		local lat1 y_cent
		
		// at least one lat/lon is a variable; others can be scalar or number
		tempvar touse
		gen `touse' = 1 `if' `in'
		qui replace `touse' = 0 if mi(`lat1',`lon1')
		qui replace `touse' = 0 if `touse' == .					
			
		tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
			cos_Sqalpha uSq A B 
			
		gen double `sin_alpha1'	= sin(rad_vertex)
		gen double `cos_alpha1' = cos(rad_vertex)
		gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
		gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
		gen double `sin_U1' = `tan_U1'*`cos_U1' 
		gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
		gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
		gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
		gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
		gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
		gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
		
		// iterations
		tempvar cont
		gen `cont' = `touse'	
		foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
			tempvar `v'
			gen double ``v'' = .
			}
		qui replace `sigma' = `dist_shift'/(`b'*`A')									
		
		local iter 0
		local more 1
		while `++iter' < `maxiter' & `more' {
			qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
			qui replace `sin_sigma' = sin(`sigma') if `cont'					
			qui replace `cos_sigma' = cos(`sigma') if `cont'					
			qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
				`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
			qui replace `sigma_old' = `sigma' if `cont'
			qui replace `sigma' = `dist_shift'/(`b'*`A') + `delta_sigma' if `cont'	

			// mark out observations that have converged
			qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
			// we are done if all observations have converged
			sum `cont', meanonly
			local more = r(max)
			}		
		
		tempvar x lambda C L lambda2

		gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
		gen double _Y_s_geo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
		gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
		gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
		gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
		gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
		gen double _X_s_geo = `lambda2'*180/c(pi)
		qui replace _Y_s_geo = _Y_s_geo*180/c(pi)
		
		drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
			`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
			`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `dist_shift' `cont' `touse'
		cap drop __0000*
		
		geodist y_cent x_cent _Y_s_geo _X_s_geo, gen(dist_s)	

		drop _X _Y
		
		rename _X_s_geo _X
		rename _Y_s_geo _Y
		
		rename _j p_scale
		replace p_scale = 0 if scale == 1
		replace p_scale = p_scale - 1 if scale > 1 & !missing(scale)	
		
		sort _ID scale aux
		} // END geo		
	
	gen orig=1 if scale==1
	replace orig=0 if orig==.
	
	tempfile scale_pareas
	save `scale_pareas'
		save "`outpoly'/scale_pareas`outname'.dta", replace
		
	* Export to temporary placebo area reference files
	sum p_scale
	local sp_min = r(min)
	local sp_max = r(max)
	
	if ("`progress'" != "no") nois _dots 0, title(Generating scale placebo areas (Total: `=`sp_max'+1'))		
	forvalues sp = `sp_min'/`sp_max' {
		if ("`progress'" != "no") nois _dots `=`sp'+1' 0
		sum scale if p_scale == `sp'
		local scale_`sp' = r(mean)
		tempfile scale_parea_`sp'
		savesome _all using "`scale_parea_`sp''" if p_scale == `sp'
		}
	
********************************************************************************		

	* Identify incidents in various areas
	if ("`progress'" != "no") {
		nois di " "					
		nois di " "	
		nois _dots 0, title(Identifying cells in scale placebo areas via `egridmatch' link (Total: `=`sp_max'+1'))		
		}
	if "`egridmatch'" == "centroid" {	
		use "`celldata'", clear
		forvalues sp = `sp_min'/`sp_max' {
			if ("`progress'" != "no") nois _dots `=`sp'+1' 0
			if "`multiarea'" == "" {
				geoinpoly y_cent x_cent using "`scale_parea_`sp''", unique 	
				rename _ID in_sparea_`sp'										
				qui replace in_sparea_`sp'=1 if !missing(in_sparea_`sp')
				qui replace in_sparea_`sp'=0 if missing(in_sparea_`sp')
				gen in_act_sparea_`sp'=(in_sparea_`sp'==1 & `etime'==1) if !missing(`etime')
				}
			if "`multiarea'" != "" {
				geoinpoly y_cent x_cent using "`scale_parea_`sp''"	
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
				replace _ID = 1 if inrange(`tid',ea_start,ea_end) & !missing(_ID)
				replace _ID = 0 if !inrange(`tid',ea_start,ea_end) | missing(_ID)
				bysort `grid' `tid': egen in_act_sparea_`sp' = max(_ID)
				drop ea_* _ID
				duplicates drop `grid' `tid' in_act_sparea_`sp', force
				}
			lab var in_act_sparea_`sp' "Centroid of (grid) cell in active operation area"
			}
		}	
		
	if "`egridmatch'" == "vertex" {		
		if "`multiarea'" == "" {
			use "`cellpoly'", clear
			forvalues sp = `sp_min'/`sp_max' {
				if ("`progress'" != "no") nois _dots `=`sp'+1' 0
						
				geoinpoly _Y _X using "`scale_parea_`sp''", unique
				rename _ID aux	
				qui replace aux=1 if !missing(aux)
				qui replace aux=0 if missing(aux)
				bysort `grid': egen in_sparea_`sp'=max(aux)
				drop aux
				}
				
			drop _X _Y
			duplicates drop

			merge 1:m `grid' using "`celldata'"
				keep if _merge==3
				drop _merge
			
			forvalues sp = `sp_min'/`sp_max' {
				gen in_act_sparea_`sp'=(in_sparea_`sp'==1 & `etime'==1) if !missing(`etime')
				lab var in_act_sparea_`sp' "Vertex of (grid) cell in active operation area"
				}
			}	
		
		if "`multiarea'" != "" {
			forvalues sp = `sp_min'/`sp_max' {
				use "`cellpoly'", clear
				if ("`progress'" != "no") nois _dots `=`sp'+1' 0
				geoinpoly _Y _X using "`scale_parea_`sp''"
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
	
				rename _ID aux	
				qui replace aux=0 if missing(aux)
		
				drop if missing(_X)
				drop _X _Y
				duplicates drop
		
				reshape wide ea_start ea_end, i(`grid') j(aux)
				
				merge 1:m `grid' using "`celldata'", keepusing(`tid')
	
				gen in_act_sparea_`sp' = 0
					lab var in_act_sparea_`sp' "Vertex of (grid) cell in active operation area"

				foreach sa of local subareas {
					replace	in_act_sparea_`sp' = 1 if inrange(`tid',ea_start`sa',ea_end`sa') ///
						& !missing(ea_start`sa') & !missing(ea_end`sa')	
					}	
				drop ea_start* ea_end* _merge
				tempfile mat_`sp'
				save `mat_`sp''
				}	
			use "`celldata'", clear
			forvalues sp = `sp_min'/`sp_max' {
				merge 1:1 `grid' `tid' using `mat_`sp''
					drop _merge
				}
			}
		}			
		
	tempfile scale_pest
	save `scale_pest'	
		save "`outres'/scale_pest`outname'.dta", replace	

	* Optional / Test mapping
	* NOTE: polygon & data can't handle temp files (messed modification)
	if ("`mapdb'" != "") & ("`mapco'" != "") {
		use "`mapdb'", clear
		spmap using "`mapco'", ///
			id(_ID) fcolor(gs12) ocolor(gs12) ///
			polygon(data("`outpoly'/scale_pareas`outname'.dta") by(orig) ocolor(black cranberry) osize(vthin)) ///
			point(data("`outpoly'/scale_pareas`outname'.dta") x(x_cent) y(y_cent) ///
			size(vsmall) fcolor(cranberry))
		graph save "`outvis'/scale_pareas`outname'.gph", replace
		graph export "`outvis'/scale_pareas`outname'.png", replace
		}	

********************************************************************************		
	*** 3.2 PLACEBO AOE ESTIMATIONS ***

	tempname memhold
	tempfile scale_pout
	cap postclose `memhold'
	postfile `memhold' scale `ts' `b_seq' `se_seq' `pv_seq' `xb_seq' `xse_seq' `xpv_seq' ///
		`diag' str30(cmd method egridmatch) ///
		using "`scale_pout'"

	use "`scale_pest'", clear
	
	if ("`progress'" != "no") nois di " "					
	if ("`progress'" != "no") nois di " "					
	if "`templacebo'" != "" {
		if ("`progress'" != "no") nois _dots 0, title(Estimation for scale placebo areas ///
		(Total: `=`sp_max'+1' areas x `n_tprange' temp. shifts = `=(`sp_max'+1)*`n_tprange''))
		}
	if "`templacebo'" == ""	{
		if ("`progress'" != "no") nois _dots 0, title(Estimation for scale placebo areas (Total: `=`sp_max'+1'))
		}
		
	local s_dot = 0
	forvalues sp = `sp_min'/`sp_max' {
		preserve
			use "`scale_pareas'", clear
			qui sum scale if p_scale == `sp', meanonly
			local scale_`sp' = r(mean)
		restore	
		
		preserve
		
		if "`templacebo'" != "" {
			xtset `grid' `tid'

			tempvar aux_treated aux_t1 aux_tT t1 tT

			bysort `grid': egen `aux_treated' = max(in_act_sparea_`sp')

			sum `tid'
			local tmin = r(min)
			local tmax = r(max)
			gen `aux_t1' = in_act_sparea_`sp' if `tid' == `tmin'
			gen `aux_tT' = in_act_sparea_`sp' if `tid' == `tmax'
			bysort `grid': egen `t1' = max(`aux_t1')
			bysort `grid': egen `tT' = max(`aux_tT')
			
			local shift_suffix ""
			forvalues s = `s_start' / `s_end' {
				if `s' < 0 {
					gen in_act_sparea_`sp'_sn_`=`s' * (-1)' = F`=`s' * (-1)'.in_act_sparea_`sp'
					local shift_suffix `shift_suffix' _sn_`=`s' * (-1)'
					if "`seriesmargin'" == "zero" & !missing(in_act_sparea_`sp') ///
						replace in_act_sparea_`sp'_sn_`=`s' * (-1)' = 0 if missing(in_act_sparea_`sp'_sn_`=`s' * (-1)')
					if "`seriesmargin'" == "preserve" & !missing(in_act_sparea_`sp') {
						replace in_act_sparea_`sp'_sn_`=`s' * (-1)' = 1 if missing(in_act_sparea_`sp'_sn_`=`s' * (-1)') & `aux_treated' == 1 & `tT' == 1
						replace in_act_sparea_`sp'_sn_`=`s' * (-1)' = 0 if missing(in_act_sparea_`sp'_sn_`=`s' * (-1)') & (`aux_treated' == 0 | `tT' == 0)
						}
					}
				if `s' > 0 {
					gen in_act_sparea_`sp'_sp_`=`s'' = L`s'.in_act_sparea_`sp'
					local shift_suffix `shift_suffix' _sp_`=`s''
					if "`seriesmargin'" == "zero" & !missing(in_act_sparea_`sp') replace in_act_sparea_`sp'_sp_`=`s'' = 0 if missing(in_act_sparea_`sp'_sp_`=`s'')
					if "`seriesmargin'" == "preserve" & !missing(in_act_sparea_`sp') {
						replace in_act_sparea_`sp'_sp_`=`s'' = 1 if missing(in_act_sparea_`sp'_sp_`=`s'') & `aux_treated' == 1 & `t1' == 1
						replace in_act_sparea_`sp'_sp_`=`s'' = 0 if missing(in_act_sparea_`sp'_sp_`=`s'') & (`aux_treated' == 0 | `t1' == 0)
						}
					}
				}
			}		
		
		foreach shift in "" `shift_suffix' {
			local s_dot = `s_dot' + 1
			if "`shift'" == "" scalar ts = 0
			if strpos("`shift'","_sn_") scalar ts = `=subinstr("`shift'","_sn_","-",1)'
			if strpos("`shift'","_sp_") scalar ts = `=subinstr("`shift'","_sp_","",1)'

			if ("`complex'" != "") {
				xtset `grid' `tid'
				
				if strpos("`complex'","TS_`evariable'") {
					bysort `grid' (`tid'): gen TS_in_act_sparea_`sp'`shift' = sum(in_act_sparea_`sp'`shift') if in_act_sparea_`sp'`shift'==1
						replace TS_in_act_sparea_`sp'`shift' = 0 if TS_in_act_sparea_`sp'`shift'==. & !missing(in_act_sparea_`sp'`shift')		
						lab var TS_in_act_sparea_`sp'`shift' "# of cumulative time periods with active AoE"	
					}
				if strpos("`complex'","TC_`evariable'") {
					bysort `grid' (`tid'): gen TC_in_act_sparea_`sp'`shift' = sum(in_act_sparea_`sp'`shift') if in_act_sparea_`sp'`shift'==1
					replace TC_in_act_sparea_`sp'`shift' = 0 if in_act_sparea_`sp'`shift' == 0
					tempvar aux
					gen `aux' = TC_in_act_sparea_`sp'`shift'
						replace `aux' = L1.`aux' if `aux' == 0 & L1.`aux' != .
						replace `aux' = L1.`aux' if `aux' == TC_in_act_sparea_`sp'`shift' & L1.`aux' != .
						replace TC_in_act_sparea_`sp'`shift' = TC_in_act_sparea_`sp'`shift' - `aux' if TC_in_act_sparea_`sp'`shift' != 0
						lab var TC_in_act_sparea_`sp'`shift' "# of consecutive time periods with active AoE"	
					}
				local est_sp = subinstr(`"`estimation'"',"`evariable'","`complex'",1)
				local est_sp = subinstr(`"`est_sp'"',"`evariable'","in_act_sparea_`sp'`shift'",.)
				}

			if ("`complex'" == "") local est_sp = subinstr(`"`estimation'"',"`evariable'","in_act_sparea_`sp'`shift'",.)
			`est_sp'		
		
			matrix b=e(b)
			matrix V=e(V)
			local indepvars : colnames b
			forvalues i = 1 / `n_b' {
				local colbase: word `i' of `indepvars' //`e(indepvars)'
				local colbase `=subinstr("`colbase'","_`sp'`shift'","",.)'
				if strpos("`colbase'","#") & strpos("`colbase'","c.") local col_int `=subinstr("`colbase'","c.","co.",.)'
				tokenize `aoe_bs_sp'
				forvalues j = 1 / `n_aoe_b' {	
					if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
						scalar b`j' = b[1,`i']
						scalar se`j' = sqrt(V[`i',`i'])
						if "`go_t'" == "1" scalar pv`j' = 2 * ttail(e(df_r),abs(b`j'/se`j'))
						if "`go_z'" == "1" scalar pv`j' = 2 * (1-normal(abs(b`j'/se`j')))
						}		
					}
				if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
					tokenize `x_bs'
					forvalues j = 1 / `n_x_b' {	
						forvalues k = 1 / `n_ovrep' { 
							if "`: word `j' of `x_bs''" == "`: word `k' of `ovreport''" {
								if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
									local x_lab_`j'_ev "`colbase'"
									scalar xb`j' = b[1,`i']
									di "Beta(`x_lab_`j'_ev') = "xb`j'
									scalar xse`j' = sqrt(V[`i',`i'])
									if "`go_t'" == "1" scalar xpv`j' = 2 * ttail(e(df_r),abs(xb`j'/xse`j'))
									if "`go_z'" == "1" scalar xpv`j' = 2 * (1-normal(abs(xb`j'/xse`j'))) //prob with logit & other var
									}					
								}
							}
						}	
					}	
				}	
			scalar scale = `scale_`sp''
			foreach j of local diag {
					if !inlist("`j'","p_F","p_chi2") scalar `j' = e(`j')
					if "`j'" == "p_F" scalar `j' = Ftail(e(df_m),e(df_r),e(F))
					if "`j'" == "p_chi2" scalar `j' = chi2tail(e(df_m),e(chi2))
					}	

			local cmd = e(cmd)
			post `memhold' (scale) `tsb' `b_seqb' `se_seqb' `pv_seqb' `xb_seqb' `xse_seqb' `xpv_seqb' ///
				`diagb' ("`cmd'") ("`method'") ("`egridmatch'")	
			
			if !inlist("`progress'","no","detail") nois _dots `s_dot' 0	
			if ("`progress'" == "detail") {
				local dot_return 0
				forvalues j = 1 / `n_aoe_b' {	
					if b`j' == 0  local dot_return 1
					}
				nois _dots `s_dot' `dot_return'
				}
			} 
		restore	
		}	
	postclose `memhold'
	if ("`progress'" != "no") nois di " "					

********************************************************************************	
	*** 3.3 PLACEBO AOE VISUALISATION ***
	
	use "`scale_pout'", clear	

	// Note: If b = 0 & se = 0 the area is too small to include any incident, and thus, 
	//			the corresponding placebo variable was ommitted from the estimation.
	if ("`evreport'" == "" | "`complex'" == "") {
		local b b1
		local se se1
		local pv pv1
		local b_title `evariable'
		}
		
	forvalues j = 1 / `n_aoe_b' {	
		gen aux = 1 if b`j' == 0 & se`j' == 0 & missing(pv`j')
		qui replace b`j' = . if aux == 1
		qui replace se`j' = . if aux == 1
		lab var b`j' "b [`aoe_lab_`j'_ev']"
		lab var se`j' "se [`aoe_lab_`j'_ev']"
		lab var pv`j' "p-value [`aoe_lab_`j'_ev']"
		drop aux	
		if ("`evreport'" != "" & "`complex'" != "") {
			if "`evreport'" == "`aoe_lab_`j'_ev'" {
				local b b`j'
				local se se`j'
				local b_title `aoe_lab_`j'_ev'
				}
			}
		}
	if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
		forvalues j = 1 / `n_x_b' {	
			cap lab var xb`j' "b [`x_lab_`j'_ev']"
			cap lab var xse`j' "se [`x_lab_`j'_ev']"
			cap lab var xpv`j' "p-value [`x_lab_`j'_ev']"
			cap format xpv`j' %9.4f
			}
		}
		
	gen `b'_cl95 = `b' - invttail(df_r,(1-95/100)/2) * `se'
	gen `b'_cu95 = `b' + invttail(df_r,(1-95/100)/2) * `se'

	sort scale	
		save "`outres'/scale_pout`outname'.dta", replace
		
	if "`templacebo'" == "" {
		graph twoway (rarea `b'_cl95 `b'_cu95 scale, color(gs14`c_damp80')) (line `b' scale, color(black)),  ///
			yline(0) xline(1) ylabel(, nogrid) xtick(0 (`scale_s') `scale_ub') xlabel(0 (0.5) `scale_ub') /// 
			xtitle("scale factor") legend(order(2 "{&beta}" 1 "95% CI")) ///
			title("AoE scale placebo estimates for {bf:`b_title'}", size(medium) color(black)) ///
			graphregion(color(white)) bgcolor(white) name(scale_gr, replace)
		}
		
	if "`templacebo'" != "" {
		local gr_stack ""
		tempvar y
		gen `y' = 0	
		levelsof ts if inrange(ts,`sr_start',`sr_end'), local(shift)
		if ("`progress'" != "no")  {
			nois di " "
			nois di " "
			nois di as text "Creating temporal AoE scale placebo diagnostic graphs"
			}
			
		if "`tpgtype'" == "line" {
			foreach s of local shift {
				if `s' < 0 {
					local sadj n
					local s1 *(-1)
					}
				if `s' > 0 {
					local sadj p
					local s1 ""
					}
				local plc white
				local ylc gs14
				if `s' == 0 {
					local plc gs14
					local ylc white
					}
				graph twoway (line `y' scale, color(cranberry) lwidth(vthin)) ///
					(rarea `b'_cl95 `b'_cu95 scale if ts == `s', color(`ylc'`c_damp80')) ///
					(line `b' scale if ts == `s', color(black)) ///
					(line `b' scale if cmd == "wtf", yaxis(2)) , ///
					xline(1, lcolor(cranberry)) ///
					ylab(, nogrid /*nolabels noticks*/ alternate) /*yscale(lwidth(none))*/ ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal)) ///
					xlab("") xscale(range(0 `scale_ub') lwidth(none)) /// 
					ytitle("") xtitle("") legend(off) ///
					graphregion(color(`plc')) plotr(color(`plc')) bgcolor(gs14) name(scale_placebo_ts_`sadj'`=`s'`s1'', replace) nodraw
				local gr_stack `gr_stack' scale_placebo_ts_`sadj'`=`s'`s1''
				}
			// fake scale label graph
			graph twoway (line `y' scale, lwidth(none)) ///
				(line `b' scale if cmd == "wtf", yaxis(2)), ///
				ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal) color(gs14)) ///
				ylabel(, nogrid /*nolabel*/ tlcolor(gs14) labcolor(gs14) alternate) yscale(lwidth(none)) ytitle("") ///
				xtick(0 (`scale_s') `scale_ub') xlabel(0 (0.5) `scale_ub') /// 
				yscale(axis(2) lwidth(none)) xtitle("scale factor") legend(off) ///
				graphregion(color(gs14)) plotr(color(gs14)) bgcolor(white) name(scale_placebo_ts_lab, replace) nodraw
				}
			
		if "`tpgtype'" != "line" {
			foreach s of local shift {
				if `s' < 0 {
					local sadj n
					local s1 *(-1)
					}
				if `s' > 0 {
					local sadj p
					local s1 ""
					}
				local plc white
				local ylc gs14
				if `s' == 0 {
					local plc gs14
					local ylc white
					}
				graph twoway (bar `b' scale if ts == `s' & `b' > 0 & `pv' < 0.01, color(`c_p1') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & `b' > 0 & `pv' < 0.05 & `pv' >= 0.01, color(`c_p5') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & `b' > 0 & `pv' < 0.1 & `pv' >= 0.05, color(`c_p10') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & `b' < 0 & `pv' < 0.01, color(`c_n1') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & `b' < 0 & `pv' < 0.05, color(`c_n5') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & `b' < 0 & `pv' < 0.1 & `pv' >= 0.05, color(`c_n10') lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(bar `b' scale if ts == `s' & (`b' == 0 | (`pv' >= 0.05 & !missing(`pv'))), color(gs12) lcolor(white) barwidth(`=0.9*`scale_s'') lwidth(none)) ///
					(line `b' scale if cmd == "wtf", yaxis(2))  ///
					(line `y' scale, color(`ylc') lwidth(thin)), ///
					/*xline(0, lcolor(`plc'))*/ xline(1, lcolor(black)) xscale(range(0 `scale_ub') lwidth(none)) ///
					ylab(, nogrid /*nolabels noticks*/ alternate) ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal)) ///
					xlab("") ytitle("") xtitle("") legend(off) ///
					graphregion(color(`plc')) plotr(color(`plc')) bgcolor(gs14) name(scale_placebo_ts_`sadj'`=`s'`s1'', replace) nodraw
				local gr_stack `gr_stack' scale_placebo_ts_`sadj'`=`s'`s1''
				}
			// fake scale label graph
			graph twoway (bar `y' scale, barwidth(`=0.9*`scale_s'') lwidth(none)) ///
				(line `b' scale if cmd == "wtf", yaxis(2)), ///
				ytitle("temporal shift:" "`s' periods", axis(2) orientation(horizontal) color(gs14)) ///
				ylabel(, nogrid tlcolor(gs14) labcolor(gs14) alternate) yscale(lwidth(none)) ytitle("") ///
				xtick(0 (`scale_s') `scale_ub') xlabel(0 (0.5) `scale_ub') /// 
				yscale(axis(2) lwidth(none)) xtitle("scale factor") legend(off) ///
				graphregion(color(gs14)) plotr(color(gs14)) bgcolor(white) name(scale_placebo_ts_lab, replace) nodraw
			}
		local n_gr: word count `gr_stack'	
		local gr_stack_in ""
		tokenize `gr_stack'
		forvalues i = 1 / `n_gr' {
			local gr_stack_in `gr_stack_in' ``=`n_gr'-`i'+1''
			}
		if ("`progress'" != "no") nois di as text "Combining temporal AoE scale placebo diagnostics graphs"
		graph combine `gr_stack_in' scale_placebo_ts_lab, graphregion(color(gs14)) ///
			title("AoE scale placebo estimates for {bf:`b_title'}", size(medium) color(black)) ///	
			cols(1) imargin(0 0 0 0) ycommon		
		}		
		
		graph save "`outvis'/scale_pdiag`outname'.gph", replace	
		graph export "`outvis'/scale_pdiag`outname'.png", replace
		}
	}
	}
	
********************************************************************************
	
qui {	
if "`design'" == "permutation" {

	global gr_stack_perm ""
		
	tempname memhold
	cap postclose `memhold'
	tempfile perm_out
	postfile `memhold' rep x_j y_j pos_radius rotation scale `b_seq' `se_seq' `pv_seq' ///
		`xb_seq' `xse_seq' `xpv_seq' `diag' str30(cmd method egridmatch seed)  ///
		using "`perm_out'"

	if ("`progress'" != "no") nois di " "		
	if ("`progress'" != "no") nois _dots 0, title(Estimation for observed and permutation placebo areas (Total: 1+`replications'))

	if "`seed'" != "" set seed `seed'
	if "`seed'" == "" local seed random
	
	local s_dot = 0
	forvalues rep = 0 / `replications' {
		local s_dot = `s_dot' + 1
		if `rep' > 0 {
			use "`area_geo'", clear	
			
			* Get random rotation point (within the predefined rotation bound)
			if "`r_distr'"=="uniform" local rnd_rot runiform(-`rot_b',`rot_b')	
			if "`r_distr'"=="centered" local rnd_rot (0.3*`rot_b')*invnorm((normal((`rot_b'-0)/(0.3*`rot_b')) ///
				-normal((-`rot_b'-0)/(0.3*`rot_b')))*runiform()+normal((-`rot_b'-0)/(0.3*`rot_b')))+0
			if `rot_b'==0 local rnd_rot 0
			gen double rotation=`rnd_rot' in 1
				replace rotation=rotation[_n-1] if missing(rotation)
			gen double rota=rotation*c(pi)/180		
			
			* Get random scaling factor (within 0, or the lower bound, and the upper bound, symmetric around 1)
			if "`s_distr'"=="uniform" local rnd_scale runiform(0,2)
			if "`s_distr'"=="centered" local rnd_scale 0.3*invnorm((normal((2-1)/0.3)-normal((0-1)/0.3))*runiform()+normal((0-1)/0.3))+1
			gen scale=`rnd_scale' in 1
				replace scale=(`scale_ub'-1)*(scale-1)+1 if scale>1
				if "`scale_lb'" != "" replace scale = scale * (1 - `scale_lb') + `scale_lb' if scale < 1	
				replace scale = scale[_n-1] if missing(scale)
				
			if "`method'" == "degree" {		
				* Get random centroid (within the jitter radius)
				if "`p_distr'"=="uniform" local rnd_jit sqrt(runiform(0,`=`pos_r'^2'))
				if "`p_distr'"=="centered" local rnd_jit  = abs((0.3*`pos_r')*invnorm((normal((`pos_r'-0)/(0.3*`pos_r')) ///
					-normal((-`pos_r'-0)/(0.3*`pos_r')))*runiform()+normal((-`pos_r'-0)/(0.3*`pos_r')))+0)	
				gen double r=`rnd_jit' in 1
					replace r=r[_n-1] if missing(r)
				gen double a=runiform(0,2*c(pi)) in 1
					replace a=a[_n-1] if missing(a)
				
				gen double pos_radius = r
				
				gen double x_j=r*cos(a)
				gen double y_j=r*sin(a)

				gen x_cent_j=x_cent+x_j	
				gen y_cent_j=y_cent+y_j	
				gen _X_j=_X+x_j	
				gen _Y_j=_Y+y_j	
				
				* Apply rotation
				gen x_jr=_X_j-x_cent_j
				gen y_jr=_Y_j-y_cent_j
				
				gen _X_jr=-y_jr*sin(-rotation*_pi/180)+x_jr*cos(-rotation*c(pi)/180)
					replace _X_jr=_X_jr+x_cent_j
				gen _Y_jr=y_jr*cos(-rotation*_pi/180)+x_jr*sin(-rotation*c(pi)/180)
					replace _Y_jr=_Y_jr+y_cent_j				
					
				* Derive vectors from centroid to vertices
				gen x_vec=_X_jr-x_cent_j
				gen y_vec=_Y_jr-y_cent_j

				// obtain rescaled vertices
				gen _X_jrs=scale*x_vec+x_cent_j
				gen _Y_jrs=scale*y_vec+y_cent_j
				} // end of degree-procedure		
				
			if "`method'" == "geodetic" {		
				local pos_r100 = `pos_r' * 100

				tempname d2r a b f
				scalar `d2r' = c(pi) / 180
				scalar `a' = 6378137
				scalar `f' = 1/298.257223563 // WGS-84 ellipsiod	
				scalar `b' = `a' - `a' * `f'
				local maxiter 1000
		
			*** JITTERING: random uniform jitter centroid selection ***
				if "`p_distr'"=="uniform" {
				
				local lon1 x_cent
				local lat1 y_cent
				
				// at least one lat/lon is a variable; others can be scalar or number
				tempvar touse
				gen `touse' = 1 `if' `in'
				qui replace `touse' = 0 if mi(`lat1',`lon1')
				qui replace `touse' = 0 if `touse' == .		
				
				// rectangular coordinate boundaries, based on four bearings (0, 90, 180, 270 degrees)
				// bound=1: north, bound=2: east, bound=3: south, bound=4: west
				// for large polygons (many vertices), more efficient calculation possible
				forvalues recbound=1/4 {
					tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
						cos_Sqalpha uSq A B 
						
					gen double `sin_alpha1'	= sin(`d2r'*90*(`recbound'-1))
					gen double `cos_alpha1' = cos(`d2r'*90*(`recbound'-1))
					gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
					gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
					gen double `sin_U1' = `tan_U1'*`cos_U1' 
					gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
					gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
					gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
					gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
					gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
					gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
					
					// iterations
					tempvar cont
					gen `cont' = `touse' in 1	
					foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
						tempvar `v'
						gen double ``v'' = .
						}
					qui replace `sigma' = `pos_r100'*1000/(`b'*`A')									
					
					local iter 0
					local more 1
					while `++iter' < `maxiter' & `more' {
						qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
						qui replace `sin_sigma' = sin(`sigma') if `cont'					
						qui replace `cos_sigma' = cos(`sigma') if `cont'					
						qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
							`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
						qui replace `sigma_old' = `sigma' if `cont'
						qui replace `sigma' = `pos_r100'*1000/(`b'*`A') + `delta_sigma' if `cont'	
			
						// mark out observations that have converged
						qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
						// we are done if all observations have converged
						sum `cont', meanonly
						local more = r(max)
						}		
					
					tempvar x lambda C L lambda2

					gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
					gen double lat_rb`recbound' = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
					gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
					gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
					gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
					gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
					gen double lon_rb`recbound' = `lambda2'*180/c(pi)
					qui replace lat_rb`recbound' = lat_rb`recbound'*180/c(pi)
					
					drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
						`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
						`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont'
					}	
					
				// required for rectangular selection: lat_rb1 lon_rb2 lat_rb3 lon_rb4
				// uniform draw
				gen x_cent_jgeo=runiform(lon_rb4,lon_rb2) in 1
				gen y_cent_jgeo=runiform(lat_rb3,lat_rb1) in 1
				
				// get distance from original centroid and geo-jittered centroid
				geodist y_cent x_cent y_cent_jgeo x_cent_jgeo, gen(geojdist)
				
				// replace out of circle coordinates until coordinates lie within
				local iter 0
				local more 1
				while `++iter' < `maxiter' & `more' {		
					replace x_cent_jgeo=runiform(lon_rb4,lon_rb2) if geojdist>`pos_r100' in 1
					replace y_cent_jgeo=runiform(lat_rb3,lat_rb1) if geojdist>`pos_r100' in 1
					drop geojdist
					geodist y_cent x_cent y_cent_jgeo x_cent_jgeo, gen(geojdist)
					sum geojdist
					if r(max)<=`pos_r100' local more = 0		
					}
				rename geojdist pos_radius	
				}	
			
		*** JITTERING: random center-clustered jitter centroid selection
			if "`p_distr'" == "centered" {
				local lon1 x_cent // to be checked
				local lat1 y_cent // to be checked
				
				tempvar touse
				gen `touse' = 1 `if' `in'
				qui replace `touse' = 0 if mi(`lat1',`lon1')
				qui replace `touse' = 0 if `touse' == .	

				gen radius_jit=runiform(0,`pos_r100')
				gen radian_jit=runiform(0,360)
				tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
					cos_Sqalpha uSq A B 
					
				gen double `sin_alpha1'	= sin(`d2r'*radian_jit)
				gen double `cos_alpha1' = cos(`d2r'*radian_jit)
				gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
				gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
				gen double `sin_U1' = `tan_U1'*`cos_U1' 
				gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
				gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
				gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
				gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
				gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
				gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
				
				// iterations
				tempvar cont
				gen `cont' = `touse' in 1	
				foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
					tempvar `v'
					gen double ``v'' = .
					}
				replace `sigma' = radius_jit*1000/(`b'*`A')	in 1	
				
				local iter 0
				local more 1
				while `++iter' < `maxiter' & `more' {
					qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
					qui replace `sin_sigma' = sin(`sigma') if `cont'					
					qui replace `cos_sigma' = cos(`sigma') if `cont'					
					qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
						`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
					qui replace `sigma_old' = `sigma' if `cont'
					qui replace `sigma' = radius_jit*1000/(`b'*`A') + `delta_sigma' if `cont'	
		
					// mark out observations that have converged
					qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
					// we are done if all observations have converged
					sum `cont', meanonly
					local more = r(max)
					}		
				
				tempvar x lambda C L lambda2

				gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
				gen double y_cent_jgeo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
				gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
				gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
				gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
				gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
				gen double x_cent_jgeo = `lambda2'*180/c(pi)
				qui replace y_cent_jgeo = y_cent_jgeo*180/c(pi)
				
				drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
					`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
					`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont'
				cap drop __0000*
				
				geodist y_cent x_cent y_cent_jgeo x_cent_jgeo, gen(pos_radius)
				}
				
			sort aux
			if "`fixpoint'" != "split" {
				replace x_cent_jgeo=x_cent_jgeo[_n-1] if x_cent_jgeo==.
				replace y_cent_jgeo=y_cent_jgeo[_n-1] if y_cent_jgeo==.
				}
				
			if "`fixpoint'" == "split" {
				replace x_cent_jgeo=x_cent_jgeo[_n-1] if x_cent_jgeo==. & _ID == _ID[_n-1]
				replace y_cent_jgeo=y_cent_jgeo[_n-1] if y_cent_jgeo==. & _ID == _ID[_n-1]
				
				* Reposition other areas according to the first areas random repositioning
					local lon1 x_cent
					local lon2 x_cent_jgeo
					local lat1 y_cent
					local lat2 y_cent_jgeo
					
					* Get radians between original centroid of one area and its random new centroid
					tempvar touse
					gen `touse' = 1 `if' `in'
					qui replace `touse' = 0 if mi(`lat1',`lon1',`lat2',`lon2')
					qui replace `touse' = 0 if `touse' == .			

					* First approximation, equation 13
					tempvar L lambda
					gen double `L' =  `d2r' * (`lon2' - `lon1') if `touse'
					gen double `lambda' = `L' if `touse'			
					
					tempvar U1 U2 sin_U1 sin_U2 cos_U1 cos_U2
					gen double `U1' = atan((1-`f') * tan(`d2r' *`lat1')) if `touse'
					gen double `U2' = atan((1-`f') * tan(`d2r' *`lat2')) if `touse'
					gen double `sin_U1' = sin(`U1') if `touse'
					gen double `sin_U2' = sin(`U2') if `touse'
					gen double `cos_U1' = cos(`U1') if `touse'
					gen double `cos_U2' = cos(`U2') if `touse'
					drop `U1' `U2'

					tempvar cont
					gen `cont' = `touse'
					foreach v in sin_sigma cos_sigma sigma sin_alpha ///
						cos_sq_alpha cos_2sigma_m C lambda_old {
						tempvar `v'
						gen double ``v'' = .
						}
						
					local iter 0
					local more 1
					while `++iter' < `maxiter' & `more' {
						// equation 14
						qui replace `sin_sigma' = sqrt((`cos_U2' * sin(`lambda'))^2 + ///
							(`cos_U1' * `sin_U2' - `sin_U1' * `cos_U2' * cos(`lambda'))^2) ///
							if `cont'
						// mark out co-incident points
						qui replace `cont' = 0 if `sin_sigma' == 0
						// equation 15
						qui replace `cos_sigma' = `sin_U1' * `sin_U2' + `cos_U1' * ///
							`cos_U2' * cos(`lambda') if `cont'
						// equation 16
						qui replace `sigma' = atan2(`sin_sigma',`cos_sigma') if `cont'
						// equation 17
						qui replace `sin_alpha' = `cos_U1' * `cos_U2' * sin(`lambda') / ///
							`sin_sigma' if `cont'
						// use trig identity to obtain cos^2 alpha
						qui replace `cos_sq_alpha' = 1 - `sin_alpha'^2 if `cont'
						// equation 18
						qui replace `cos_2sigma_m' = `cos_sigma' - 2 * `sin_U1' * ///
							`sin_U2' / `cos_sq_alpha' if `cont'
						// adjust if both points are on the equator
						qui replace `cos_2sigma_m' = 0 if `cos_sq_alpha' == 0 & `cont'
						// compute new lambda and compare to previous one
						qui replace `lambda_old' = `lambda' if `cont'
						// equation 10
						qui replace `C' = `f' / 16 * `cos_sq_alpha' * ///
							(4 + `f' * (4 - 3 * `cos_sq_alpha')) if `cont'
						// equation 11
						qui replace `lambda' = `L' + (1 - `C') * `f' * `sin_alpha' * ///
						  (	`sigma' + `C' * `sin_sigma' * (`cos_2sigma_m' + ///
						  `C'*`cos_sigma' * (-1 + 2* `cos_2sigma_m'^2))) if `cont'
						// mark out observations that have converged
						qui replace `cont' = 0 if abs(`lambda'-`lambda_old') <= 1e-12
						// we are done if all observations have converged
						sum `cont', meanonly
						local more = r(max)
						}
					
						// Radians between original centroid and first random centroid
						gen double rad_1repos=atan2(`cos_U2'*sin(`lambda'),`cos_U1'*`sin_U2'-`sin_U1'*`cos_U2'*cos(`lambda'))
							*replace rad_1repos=mod(rad_1repos+2*c(pi), 2*c(pi))
							
						replace rad_1repos = rad_1repos[_n-1] if missing(rad_1repos)			
						replace pos_radius = pos_radius[_n-1] if missing(pos_radius)			
						
						* Shift other areas' centroids
						local lon1 x_cent
						local lat1 y_cent

						// at least one lat/lon is a variable; others can be scalar or number
						tempvar touse
						gen `touse' = 1 `if' `in'
						qui replace `touse' = 0 if mi(`lat1',`lon1')
						qui replace `touse' = 0 if `touse' == .					
						
						tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
							cos_Sqalpha uSq A B 
							
						gen double `sin_alpha1'	= sin(rad_1repos)
						gen double `cos_alpha1' = cos(rad_1repos)
						gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
						gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
						gen double `sin_U1' = `tan_U1'*`cos_U1' 
						gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
						gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
						gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
						gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
						gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
						gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
						
						// iterations
						tempvar cont
						gen `cont' = `touse'	
						foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
							tempvar `v'
							gen double ``v'' = .
							}
						qui replace `sigma' = pos_radius*1000/(`b'*`A')									
						
						local iter 0
						local more 1
						while `++iter' < `maxiter' & `more' {
							qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
							qui replace `sin_sigma' = sin(`sigma') if `cont'					
							qui replace `cos_sigma' = cos(`sigma') if `cont'					
							qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
								`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
							qui replace `sigma_old' = `sigma' if `cont'
							qui replace `sigma' = pos_radius*1000/(`b'*`A') + `delta_sigma' if `cont'	
				
							// mark out observations that have converged
							qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
							// we are done if all observations have converged
							sum `cont', meanonly
							local more = r(max)
							}		
					
						tempvar x lambda C L lambda2 
						gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
						replace y_cent_jgeo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))*180/c(pi) ///
							if missing(y_cent_jgeo)
						gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
						gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
						gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
						gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
						replace x_cent_jgeo = `lambda2'*180/c(pi) if missing(x_cent_jgeo)
						
						drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
							`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
							`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont' `touse'
						cap drop __0000*	
						}
				
			// change in coordinate degrees (for reporting only)
			gen x_j= x_cent_jgeo - x_cent
			gen y_j= y_cent_jgeo - y_cent
		
			*** JITTERING: Get new vertices in relation to jittered centroid ***
			
			local lon1 x_cent_jgeo
			local lat1 y_cent_jgeo

			// at least one lat/lon is a variable; others can be scalar or number
			tempvar touse
			gen `touse' = 1 `if' `in'
			qui replace `touse' = 0 if mi(`lat1',`lon1')
			qui replace `touse' = 0 if `touse' == .					
			
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(rad_vertex)
			gen double `cos_alpha1' = cos(rad_vertex)
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			// iterations
			tempvar cont
			gen `cont' = `touse'	
			foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
				tempvar `v'
				gen double ``v'' = .
				}
			qui replace `sigma' = disto*1000/(`b'*`A')									
			
			local iter 0
			local more 1
			while `++iter' < `maxiter' & `more' {
				qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
				qui replace `sin_sigma' = sin(`sigma') if `cont'					
				qui replace `cos_sigma' = cos(`sigma') if `cont'					
				qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
					`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
				qui replace `sigma_old' = `sigma' if `cont'
				qui replace `sigma' = disto*1000/(`b'*`A') + `delta_sigma' if `cont'	

				// mark out observations that have converged
				qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
				// we are done if all observations have converged
				sum `cont', meanonly
				local more = r(max)
				}		
			
			tempvar x lambda C L lambda2

			gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
			gen double _Y_j_geo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
			gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
			gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
			gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
			gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
			gen double _X_j_geo = `lambda2'*180/c(pi)
			qui replace _Y_j_geo = _Y_j_geo*180/c(pi)
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
				`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont' `touse'
			cap drop __0000*

			// vertex coordinates not properly shifted
			geodist y_cent_jgeo x_cent_jgeo _Y_j_geo _X_j_geo, gen(dist_j)				
			
		*** ROTATING: Get new vertices in relation to jittered centroid & random rotation ***

			// Add random rotation component to original radian
			gen double rnd_rad_vertex=rad_vertex+rota
		
			local lon1 x_cent_jgeo
			local lat1 y_cent_jgeo

			// at least one lat/lon is a variable; others can be scalar or number
			tempvar touse
			gen `touse' = 1 `if' `in'
			qui replace `touse' = 0 if mi(`lat1',`lon1')
			qui replace `touse' = 0 if `touse' == .					
			
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(rnd_rad_vertex)
			gen double `cos_alpha1' = cos(rnd_rad_vertex)
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			// iterations
			tempvar cont
			gen `cont' = `touse'	
			foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
				tempvar `v'
				gen double ``v'' = .
				}
			qui replace `sigma' = disto*1000/(`b'*`A')									
			
			local iter 0
			local more 1
			while `++iter' < `maxiter' & `more' {
				qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
				qui replace `sin_sigma' = sin(`sigma') if `cont'					
				qui replace `cos_sigma' = cos(`sigma') if `cont'					
				qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
					`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
				qui replace `sigma_old' = `sigma' if `cont'
				qui replace `sigma' = disto*1000/(`b'*`A') + `delta_sigma' if `cont'	

				// mark out observations that have converged
				qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
				// we are done if all observations have converged
				sum `cont', meanonly
				local more = r(max)
				}		
			
			tempvar x lambda C L lambda2

			gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
			gen double _Y_jr_geo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
			gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
			gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
			gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
			gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
			gen double _X_jr_geo = `lambda2'*180/c(pi)
			qui replace _Y_jr_geo = _Y_jr_geo*180/c(pi)
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
				`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `cont' `touse'
			cap drop __0000*

			// vertex coordinates not properly shifted
			geodist y_cent_jgeo x_cent_jgeo _Y_jr_geo _X_jr_geo, gen(dist_jr)			
			
		*** SCALING: get vertices for scaled (jittered & rotated) polygon ***

			tempvar dist_shift
			gen double `dist_shift' = 1000*disto*scale

			local lon1 x_cent_jgeo
			local lat1 y_cent_jgeo
			
			// at least one lat/lon is a variable; others can be scalar or number
			tempvar touse
			gen `touse' = 1 `if' `in'
			qui replace `touse' = 0 if mi(`lat1',`lon1')
			qui replace `touse' = 0 if `touse' == .					
				
			tempvar sin_alpha1 cos_alpha1 tan_U1 cos_U1 sin_U1 sigma1 sin_alpha ///
				cos_Sqalpha uSq A B 
				
			gen double `sin_alpha1'	= sin(rnd_rad_vertex)
			gen double `cos_alpha1' = cos(rnd_rad_vertex)
			gen double `tan_U1' = (1-`f')*tan(`d2r'*`lat1') 		
			gen double `cos_U1' = 1/sqrt(1+`tan_U1'^2) 
			gen double `sin_U1' = `tan_U1'*`cos_U1' 
			gen double `sigma1' = atan2(`tan_U1',`cos_alpha1') 		
			gen double `sin_alpha' = `cos_U1'*`sin_alpha1' 			
			gen double `cos_Sqalpha' = 1 - `sin_alpha'^2
			gen double `uSq' = `cos_Sqalpha'*(`a'^2-`b'^2)/(`b'^2)	
			gen double `A' = 1+`uSq'/16384*(4096+`uSq'*(-768+`uSq'*(320-175*`uSq')))	
			gen double `B' = `uSq'/1024*(256+`uSq'*(-128+`uSq'*(74-47*`uSq')))			
			
			// iterations
			tempvar cont
			gen `cont' = `touse'	
			foreach v in cos_2sigmaM sin_sigma cos_sigma delta_sigma sigma sigma_old  {
				tempvar `v'
				gen double ``v'' = .
				}
			qui replace `sigma' = `dist_shift'/(`b'*`A')									
			
			local iter 0
			local more 1
			while `++iter' < `maxiter' & `more' {
				qui replace `cos_2sigmaM' = cos(2*`sigma1'+`sigma') if `cont'		
				qui replace `sin_sigma' = sin(`sigma') if `cont'					
				qui replace `cos_sigma' = cos(`sigma') if `cont'					
				qui replace `delta_sigma' = `B'*`sin_sigma'*(`cos_2sigmaM'+`B'/4*(`cos_sigma'*(-1+2*`cos_2sigmaM'^2)- ///
					`B'/6*`cos_2sigmaM'*(-3+4*`sin_sigma'^2)*(-3+4*`cos_2sigmaM'^2))) if `cont'			
				qui replace `sigma_old' = `sigma' if `cont'
				qui replace `sigma' = `dist_shift'/(`b'*`A') + `delta_sigma' if `cont'	

				// mark out observations that have converged
				qui replace `cont' = 0 if abs(`sigma'-`sigma_old') <= 1e-12
				// we are done if all observations have converged
				sum `cont', meanonly
				local more = r(max)
				}		
			
			tempvar x lambda C L lambda2

			gen double `x' = `sin_U1'*`sin_sigma'-`cos_U1'*`cos_sigma'*`cos_alpha1'	
			gen double _Y_jrs_geo = atan2(`sin_U1'*`cos_sigma'+`cos_U1'*`sin_sigma'*`cos_alpha1',(1-`f')*sqrt(`sin_alpha'^2+`x'^2))	
			gen double `lambda' = atan2(`sin_sigma'*`sin_alpha1',`cos_U1'*`cos_sigma'-`sin_U1'*`sin_sigma'*`cos_alpha1') 	
			gen double `C' = `f'/16*`cos_Sqalpha'*(4+`f'*(4-3*`cos_Sqalpha'))		
			gen double `L' = `lambda'-(1-`C')*`f'*`sin_alpha'*(`sigma'+`C'*`sin_sigma'*(`cos_2sigmaM'+`C'*`cos_sigma'*(-1+2*`cos_2sigmaM'^2))) 
			gen double `lambda2' =  mod(`d2r'*`lon1'+`L'+3*c(pi),2*c(pi))-c(pi)
			gen double _X_jrs_geo = `lambda2'*180/c(pi)
			qui replace _Y_jrs_geo = _Y_jrs_geo*180/c(pi)
			
			drop `sin_alpha1' `cos_alpha1' `tan_U1' `cos_U1' `sin_U1' `sigma1' `sin_alpha' ///
				`cos_Sqalpha' `uSq' `A' `B' `x' `lambda' `C' `L' `lambda2' `cos_2sigmaM' ///
				`sin_sigma' `cos_sigma' `delta_sigma' `sigma' `sigma_old' `dist_shift' `cont' `touse'
			cap drop __0000*
			
			geodist y_cent_jgeo x_cent_jgeo _Y_jrs_geo _X_jrs_geo, gen(dist_jrs)	

			} // end of geo-distance procedure

			// Streamline resulting data set
			rename _X _X_0
			rename _Y _Y_0
			
			reshape long _X _Y, i(_ID aux) string
				replace _j = subinstr(_j,"_geo","",.)

			gen type=0 if _j=="_0"
			replace type=1 if _j=="_j"
			replace type=2 if _j=="_jr"
			replace type=3 if _j=="_jrs"
			lab def type 0 "original" 1 "re-positioned" 2 "re-positioned & rotated" ///
				3 "re-positioned & rotated & re-scaled", modify
			lab val type type	
			
			drop _j
			
			keep _ID aux x_cent y_cent rota scale x_j y_j pos_radius _Y _X type		
			
			replace rota=rota*180/c(pi) // in degrees
			rename rota rotation
			
			sort _ID type aux
			
			if inrange(`rep',1,6) {	
				replace rotation = 0 if inrange(type,0,1)
				replace scale = 1 if inrange(type,0,2)
				replace x_j = 0 if type == 0
				replace y_j = 0 if type == 0
				replace x_cent = x_cent + x_j
				replace y_cent = y_cent + y_j
				save "`outpoly'/perm_areas_rep`rep'`outname'.dta", replace
				}
			
			keep if type == 3
			foreach d in x_j y_j pos_radius rotation scale {
				sum `d', meanonly
				local `d' = r(mean) 
				local `d'_di `: di %9.2f r(mean)'
				}
				
			keep _ID _Y _X type // tmp
			tempfile perm_areas
			
			save `perm_areas'
			}
			
********************************************************************************
	
	* Optional / Test mapping
	* NOTE: polygon & data can't handle temp files (messed mdoification)	
	if ("`mapdb'" != "") & ("`mapco'" != "") {
		if inrange(`rep',1,6) {	
			use "`mapdb'", clear
			spmap using "`mapco'", ///
				id(_ID) fcolor(gs12) ocolor(gs12) ///
				polygon(data("`outpoly'/perm_areas_rep`rep'`outname'.dta") by(type) ocolor(black midgreen orange red) ///
				osize(vthin) legenda(on)) legend(symysize(*0.75) rows(1) region(fcolor(white))) ///
				point(data("`outpoly'/perm_areas_rep`rep'`outname'.dta") by(type) xcoord(x_cent) ycoord(y_cent) ///
				select(keep if inlist(type,0,3)) size(tiny tiny) fcolor(black red)) ///
				title("Position: {&Delta}X=`x_j_di', {&Delta}Y=`y_j_di'" ///
				"Rotation: `rotation_di' degrees" "Scaling factor: `scale_di'", color(black) size(medium)) ///
				name(rep_`rep', replace) nodraw
			graph save "`outvis'/rep_`rep'.gph", replace		//tmp
			global gr_stack_perm $gr_stack_perm rep_`rep'
			}	
		}	
			
	* Identify incidents in various areas
		if "`egridmatch'" == "centroid" {	
			use using "`celldata'", clear		

			if "`multiarea'" == "" {
				if `rep' == 0 {
					geoinpoly y_cent x_cent using "`areadata'", unique	
					}
				if `rep' > 0 {	
					geoinpoly y_cent x_cent using "`perm_areas'", unique 	
					}
				rename _ID in_parea										
				qui replace in_parea = 1 if !missing(in_parea)
				qui replace in_parea = 0 if missing(in_parea)
				gen in_act_parea = (in_parea==1 & `etime'==1) if !missing(`etime')				
				}
				
			if "`multiarea'" != "" {
				if `rep' == 0 {
					geoinpoly y_cent x_cent using "`areadata'"	
					}
				if `rep' > 0 {	
					geoinpoly y_cent x_cent using "`perm_areas'" 	
					}
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}
				replace _ID = 1 if inrange(`tid',ea_start,ea_end) & !missing(_ID)
				replace _ID = 0 if !inrange(`tid',ea_start,ea_end) | missing(_ID)
				bysort `grid' `tid': egen in_act_parea = max(_ID)
				drop ea_* _ID
				duplicates drop `grid' `tid' in_act_parea, force		
				}				
			lab var in_act_parea "Centroid of (grid) cell in active operation area"		
			}
			
		if "`egridmatch'" == "vertex" {		
			if "`multiarea'" == "" {
				use "`cellpoly'", clear			
				
				if `rep' == 0 {
					geoinpoly _Y _X using "`areadata'", unique 	
					}
				if `rep' > 0 {	
					geoinpoly _Y _X using "`perm_areas'", unique 	
					}		
					
				rename _ID aux	
				qui replace aux=1 if !missing(aux)
				qui replace aux=0 if missing(aux)
				bysort `grid': egen in_parea=max(aux)
				drop aux
				
				drop _X _Y
				duplicates drop

				merge 1:m `grid' using "`celldata'"
					keep if _merge==3
					drop _merge
				
				gen in_act_parea = (in_parea==1 & `etime'==1) if !missing(`etime')			
				}
				
			if "`multiarea'" != "" {
				use "`cellpoly'", clear
				
				if `rep' == 0 {
					geoinpoly _Y _X using "`areadata'" 	
					}
				if `rep' > 0 {	
					geoinpoly _Y _X using "`perm_areas'" 	
					}		
								
				gen ea_start = .
				gen ea_end = .
				foreach sa of local subareas {
					replace ea_start = `ea_start_`sa'' if _ID == `sa'
					replace ea_end = `ea_end_`sa'' if _ID == `sa'
					}

				rename _ID aux	
				replace aux=0 if missing(aux)
		
				drop if missing(_X)
				drop _X _Y
				duplicates drop
		
				reshape wide ea_start ea_end, i(`grid') j(aux)
				
				merge 1:m `grid' using "`celldata'" /*, keepusing(`tid')*/

				gen in_act_parea = 0

				foreach sa of local subareas {
					cap confirm var ea_start`sa' ea_end`sa'
					if !_rc {
						replace	in_act_parea = 1 if inrange(`tid',ea_start`sa',ea_end`sa') ///
							& !missing(ea_start`sa') & !missing(ea_end`sa')	
						}	
					}	
				drop ea_start* ea_end* _merge
				}	
			lab var in_act_parea "Vertex of (grid) cell in active operation area"			
			}	
		
********************************************************************************
* Estimation			
		if ("`complex'" != "") {
			xtset `grid' `tid'
				
			if strpos("`complex'","TS_`evariable'") {
				bysort `grid' (`tid'): gen TS_in_act_parea = sum(in_act_parea) if in_act_parea==1
					replace TS_in_act_parea = 0 if TS_in_act_parea==. & !missing(in_act_parea)		
					lab var TS_in_act_parea "# of cumulative time periods with active AoE"	
				}
			if strpos("`complex'","TC_`evariable'") {
				bysort `grid' (`tid'): gen TC_in_act_parea = sum(in_act_parea) if in_act_parea==1
				replace TC_in_act_parea = 0 if in_act_parea == 0
				tempvar aux
				gen `aux' = TC_in_act_parea
					replace `aux' = L1.`aux' if `aux' == 0 & L1.`aux' != .
					replace `aux' = L1.`aux' if `aux' == TC_in_act_parea & L1.`aux' != .
					replace TC_in_act_parea = TC_in_act_parea - `aux' if TC_in_act_parea != 0
					lab var TC_in_act_parea "# of consecutive time periods with active AoE"	
				}
			local est_p = subinstr(`"`estimation'"',"`evariable'","`complex'",1)
			local est_p = subinstr(`"`est_p'"',"`evariable'","in_act_parea",.)
			}	
				
			if ("`complex'" == "") local est_p = subinstr(`"`estimation'"',"`evariable'","in_act_parea",.)
			`est_p'			
			
			matrix b=e(b)
			matrix V=e(V)
			local indepvars : colnames b
			forvalues i = 1 / `n_b' {
				local colbase: word `i' of `indepvars' // `e(indepvars)'
				if strpos("`colbase'","#") & strpos("`colbase'","c.") local col_int `=subinstr("`colbase'","c.","co.",.)'
				tokenize `aoe_bs'
				forvalues j = 1 / `n_aoe_b' {
					if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
						scalar b`j' = b[1,`i']
						scalar se`j' = sqrt(V[`i',`i'])
						if "`go_t'" == "1" scalar pv`j' = 2 * ttail(e(df_r),abs(b`j'/se`j'))
						if "`go_z'" == "1" scalar pv`j' = 2 * (1-normal(abs(b`j'/se`j')))
						}		
					}
				if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
					tokenize `x_bs'
					forvalues j = 1 / `n_x_b' {	
						forvalues k = 1 / `n_ovrep' { 
							if "`: word `j' of `x_bs''" == "`: word `k' of `ovreport''" {
								if inlist("`colbase'","``j''","o``j''","o.``j''","`col_int'") {
									local x_lab_`j'_ev "`colbase'"
									scalar xb`j' = b[1,`i']
									di "Beta(`x_lab_`j'_ev') = "xb`j'
									scalar xse`j' = sqrt(V[`i',`i'])
									if "`go_t'" == "1" scalar xpv`j' = 2 * ttail(e(df_r),abs(xb`j'/xse`j'))
									if "`go_z'" == "1" scalar xpv`j' = 2 * (1-normal(abs(xb`j'/xse`j'))) //prob with logit & other var
									}					
								}
							}
						}	
					}	
				}	
			foreach j of local diag {
				if !inlist("`j'","p_F","p_chi2") scalar `j' = e(`j')
				if "`j'" == "p_F" scalar `j' = Ftail(e(df_m),e(df_r),e(F))
				if "`j'" == "p_chi2" scalar `j' = chi2tail(e(df_m),e(chi2))
				}
			if `rep' == 0 {
				local x_j = 0
				local y_j = 0
				local pos_radius = 0
				local rotation = 0			
				local scale = 1 
				}
			scalar rep = `rep'
			scalar x_j = `x_j'
			scalar y_j = `y_j'
			scalar pos_radius = `pos_radius'
			scalar rotation = `rotation'			
			scalar scale =`scale' 
			local cmd = e(cmd)
			post `memhold' (rep) (x_j) (y_j) (pos_radius) (rotation) (scale)  `b_seqb' `se_seqb' `pv_seqb' ///
				`xb_seqb' `xse_seqb' `xpv_seqb' ///
				`diagb' ("`cmd'") ("`method'") ("`egridmatch'") ("`seed'") 
				
			if !inlist("`progress'","no","detail") nois _dots `s_dot' 0	
			if ("`progress'" == "detail") {
				local dot_return 0
				forvalues j = 1 / `n_aoe_b' {	
					if b`j' == 0  local dot_return 1
					}
				nois _dots `s_dot' `dot_return'
				}
			} // REPS
		postclose `memhold'	

********************************************************************************
* Permutation visualisation	(maps)
	if ("`mapdb'" != "") & ("`mapco'" != "") {
		grc1leg $gr_stack_perm , graphregion(color(white)) cols(3)
		graph save "`outvis'/perm_pareas6`outname'.gph", replace	
		graph export "`outvis'/perm_pareas6`outname'.png", replace		
		}
		
	* Simple output visualissation	
	use "`perm_out'", clear
	
	if ("`evreport'" == "" | "`complex'" == "") {
		local b b1
		local se se1
		local pv pv1
		local b_title `evariable'
		}
		
	forvalues j = 1 / `n_aoe_b' {	
		gen aux = 1 if b`j' == 0 & se`j' == 0 & missing(pv`j')
		qui replace b`j' = . if aux == 1
		qui replace se`j' = . if aux == 1
		lab var b`j' "b [`aoe_lab_`j'_ev']"
		lab var se`j' "se [`aoe_lab_`j'_ev']"
		lab var pv`j' "p-value [`aoe_lab_`j'_ev']"
		drop aux	
		if ("`evreport'" != "" & "`complex'" != "") {
			if "`evreport'" == "`aoe_lab_`j'_ev'" {
				local b b`j'
				local se se`j'
				local pv pv`j'
				local b_title `aoe_lab_`j'_ev'
				}
			}
		}
	if "`ovreport'" != "" & `n_b' > `n_aoe_b' {
		forvalues j = 1 / `n_x_b' {	
			cap lab var xb`j' "b [`x_lab_`j'_ev']"
			cap lab var xse`j' "se [`x_lab_`j'_ev']"
			cap lab var xpv`j' "p-value [`x_lab_`j'_ev']"
			cap format xpv`j' %9.4f
			}
		}	
	lab def rep 0 "original area"
	lab val rep rep
	save "`outres'/perm_out`outname'.dta", replace
	
	if "`permreport'" != "no"  { 
	
		if "`p_bounds'" != "" {
			local pcond inrange(abs(x_j),`pbl',`pbu') & inrange(abs(y_j),`pbl',`pbu')
			local pos_res `""      Position:       |{&Delta}X| {&element} [`pbl',`pbu'], |{&Delta}Y| {&element} [`pbl',`pbu']""'
			}
			
		if "`r_bounds'" != "" {
			local rcond inrange(abs(rotation),`rbl',`rbu')
			local rot_res `""      Orientation:  |{&theta}| {&element} [`rbl',`rbu']""'
			}	
			
		if "`s_bounds'" != "" {
			local scond inrange(scale,`sbl',`sbu')
			local scale_res `""      Scale:            s {&element} [`sbl',`sbu']""'
			}

		if "`pcond'" != "" | "`rcond'" != "" | "`scond'" != "" local ifpermcond if	
		foreach cond in "`pcond'" "`rcond'" "`scond'" {
			if "`cond'" != "" local permcond `permcond' & `cond'
			if "`cond'" != "" & "`ifpermcond'" != "if" local ifpermcond `ifpermcond' & `cond'
			if "`cond'" != "" & "`ifpermcond'" == "if" local ifpermcond `ifpermcond' `cond'
			}
			
		if "`pcond'" == "" & "`rcond'" == "" & "`scond'" == "" local restr "none"	

		tempvar tz
		gen `tz' = `b'/`se'
		
		global gr_stack ""	
		foreach stat in `b' `tz' {
			count if `stat'!=. `permcond'
			tempvar y_`stat' x_`stat'
			local label_`stat': var label `stat'
			qui sum `stat' if rep == 0
			local real_`stat' = round(`r(max)',.0001)
			if `real_`stat'' < 0 {
				count if `stat' < `real_`stat'' & rep!=0 `permcond'
				local next = `r(N)'
				count if !missing(`stat') & rep!=0 `permcond'
				local nest = `r(N)'
				local plab P({&beta}{subscript:placebo} < {&beta}{subscript:real}) = `: di %9.4f `next'/`nest'' // 0`=round(`next'/`nest',.0001)'
				local tzlab P(test stat.{subscript:placebo} < test stat.{subscript:real}) = `: di %9.4f `next'/`nest'' //0`=round(`next'/`nest',.0001)'
				}
			if `real_`stat'' > 0 {
				count if `stat' > `real_`stat'' & rep!=0 `permcond'
				local next = `r(N)'
				count if !missing(`stat') & rep!=0 `permcond'
				local nest = `r(N)'
				local plab P({&beta}{subscript:placebo} > {&beta}{subscript:real}) =`:di %9.4f `next'/`nest'' //0`=round(`next'/`nest',.0001)'
				local tzlab P(test stat.{subscript:placebo} > test stat.{subscript:real}) =`:di %9.4f `next'/`nest'' //0`=round(`next'/`nest',.0001)'
				}
			_kdens `stat' `ifpermcond', gen(`y_`stat'' `x_`stat'')
			sum `y_`stat''
			local y_text = 0.5 * `r(max)'
			sum `stat' if rep!=0 `permcond'
			local reps=`r(N)'
			if `r(min)' < 0 local min = r(min)*1.001
			if `r(max)' < 0 local max = r(max)*0.999
			if `r(min)' > 0 local min = r(min)*0.999
			if `r(max)' > 0 local max = r(max)*1.001
			sum `stat' `ifpermcond'
			if `r(min)' < 0 local min_all = r(min)*1.001
			if `r(max)' < 0 local max_all = r(max)*0.999
			if `r(min)' > 0 local min_all = r(min)*0.999
			if `r(max)' > 0 local max_all = r(max)*1.001	
			local reppos = (`max_all' - `min_all')/2
			if `real_`stat'' < `=(`max_all' - `min_all')/2 + `min_all'' {
				local place ne
				if `real_`stat'' < 0 local bgap `=`real_`stat'' * 0.9'
				if `real_`stat'' > 0 local bgap `=`real_`stat'' * 1.1'
				}
			if `real_`stat'' > `=(`max_all' - `min_all')/2 + `min_all'' {
				local place nw
				if `real_`stat'' < 0 local bgap `=`real_`stat'' * 1.1'
				if `real_`stat'' > 0 local bgap `=`real_`stat'' * 0.9'			
				}
			
			if "`stat'"=="`b'" local addtxt "{&beta} =`:di %9.4f `real_`stat'''" " " "`plab'" //" " "No. of estimable coefficients: `nest'"
			if "`stat'"=="`tz'" local addtxt "test statistic ({&beta}/se) =`:di %9.4f `real_`tz'''" " " "`tzlab'" //" " "No. of estimable coefficients: `nest'"	
			
			local xlab `=round(`min',.0001)' `=round(`max',.0001)'
			if `min'< 0 & `max'> 0 local xlab `=round(`min',.0001)' 0 `=round(`max',.0001)'
		
		if ("`progress'" != "no") & "`stat'" == "`b'" {
			nois di " "
			nois di " "
			nois di as text "Drawing kernel density plots for: `b_title'"
			}
		kdens `stat' if rep!=0 `permcond', range(`min_all' `max_all') ll(`min') ul(`max') graphregion(color(white)) bgcolor(white) ///
			xline(`real_`stat'', lcolor(cranberry)) xscale(range(`min_all' `max_all')) /*xtick(`min_all' `max_all')*/ ///
			xlab(`xlab') xline(`min', lcolor(gs12) lpattern(shortdash)) xline(`max', lcolor(gs12) lpattern(shortdash)) ///
			ylabel("", nogrid) yscale(lwidth(none)) ytitle("") xtitle("") text(`y_text' `bgap' "`addtxt'", size(medlarge) orientation(horizontal) ///
			just(left) margin(small) color(cranberry) place(`place') box bcolor(white%70)) ///
			name(gr`stat', replace) nodraw 
		global gr_stack $gr_stack gr`stat'	
		}	

		graph combine $gr_stack, rows(2) graphregion(color(white) lcolor(white))  ///
			title("Overview of permutation-based coefficients" "and test statistics for {bf:`b_title'}", color(black)) ///
			note("{bf:Restrictions}:    `restr'                                                            {bf:Estimable coefficients}: `nest'" ///
			`pos_res' `rot_res' `scale_res', size(small)) iscale(0.6)
				
			graph save "`outvis'/perm_kdens`outname'.gph", replace	
			graph export "`outvis'/perm_kdens`outname'.png", replace
		}		
	} 
	}	
	
********************************************************************************
*** END OF PROGRAMM OPERATIONS & DISPLAY									 ***	
********************************************************************************	

if ("`progress'" != "no") {
	di " "
	di "*********************************************************************************"
	}
di " "	
di as text `"Placebo test visualisation graphs can be found {browse `"`outvis'"' :here}."'
*di "     `outvis'"
di " "	
di as text `"Placebo test result data can be found {browse `"`outres'"' :here}."'
*di "     `outres'"	
di " "	
di as text `"Placebo area data can be found {browse `"`outpoly'"' :here}."'
*di "     `outpoly'"	
di " "
di in red "*********************************************************************************"
di " "

********************************************************************************

end
